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MONTE CARLO CODE FOR SOLUTION OF PLANAR ELECTRON- DIODE PROBLEMS 
INCLUDING ELECTRON-NEUTRAL ELASTIC COLLISIONS 
by Paul Swigert and Charles M. Goldstein 
Lewis Research Center 

SUMMARY 

A Monte Carlo electron transport code for the self-consistent potential solutions of 
one -dimensional planar electron-diode problems including electron-neutral elastic colli- 
sions capable of employing differential scattering cross sections is presented. An ana- 
lytical description of a class of problems for which this code was written and the methods 
and techniques for solution of these problems are also presented. The code is given in- 
cluding instructions for the user, flow charts, and listings of all FORTRAN IV programs. 
Also included is material on curve and surface fits, convergence of a second-order differ- 
ential equation, machine language routines, and a sample problem. 


INTRODUCTION 

The solution of steady-state electron transport problems in the presence of a low 
density, scattering background gas and a nonuniform electric field has not as yet been 
achieved by the usual analytical methods. In the field of neutron transport problems, 
where the usual methods of analysis also fail, resort is made to the Monte Carlo method 
to obtain particular solutions. This report describes a Monte Carlo code ENEC (electron- 
neutral elastic collisions) for the self-consistent potential field solution of a class of 
electron-diode problems including the effects of electron-neutral elastic collisions. 

A preliminary version of this code, restricted to hard-sphere electron-neutral colli- 
sions was published in appendix 1 of reference 1. The present report, however, is com- 
pletely self-contained and presents many improvements. In addition, the restriction to 
hard-sphere collisions has been removed in order to treat energy- and angle -dependent 
differential scattering cross sections. 

An analytical description of the class of problems for which this code was written is 
given, and then the code itself, ENEC, is presented. The mathematical symbols used in 



the analysis are defined in appendix A. One -dimensional spline curve fits and two- 
dimensional spline surface fits are discussed in appendix B. A convergence experiment 
is described, the results obtained are discussed, and some conclusions are drawn in 
appendix C. Presented in appendix D is an explanation of an improved square root rou- 
tine. Appendix E contains machine language routines used in ENEC, which were not pro- 
gramed by the authors, but are available in the Lewis 7094 Library. The output of a 
sample problem is given in appendix F. 


ANALYSIS 

DESCRIPTION OF PROBLEM 

The geometric configuration of one -dimensional field, flux, and electrodes is depic- 
ted in figure 1. The one -dimensional problem is treated wherein the emitter and collector 



^Collector 


Figure 1. - Configuration of one-dlmenslonal field, flux, and 
electrodes. 


are assumed to be infinite parallel planes. The electric field and x-direction are normal 
to the electrode surfaces. 

The class of problems treated here is considerably more complicated than the neu- 
tron transport problems because of the nonlinearity introduced by the potential field, the 
existences of curvilinear rather than rectilinear trajectories, and the spatial as well as 
energy variation of the mean free path. 

The potential distribution is obtained as a solution to Poisson’s equation 


2 


( 1 ) 


by Picard iteration. (All symbols are defined in appendix A.) For every assumed or 
computed potential distribution ^^(x), the electron charge density p(x) is obtained by a 
Monte Carlo calculation. The question of convergence is discussed in the section 
CONVERGENCE OF SOLUTION and in appendix C. 

The code assumes that the electrons are thermionically emitted with a half- 
Maxwellian velocity distribution: 


f(u,v,w)du dv dw = 2 


m 


^TTkT, 


.3/2 -/2kTe 


0 < U < 00 


< V < oo y 


( 2 ) 


- 00 < w< 00 


Only minor modifications are necessary, however, to treat monoenergetic beam emission 
(see ref. 1). 


SOLUTION OF DIFFERENTIAL EQUATION 


In dimensionless variables, equation (1) becomes 


= c • n(x) 


dx^ 


(3) 


where n(x) is the dimensionless electron density distribution, ^(x) is the dimensionless 
potential distribution, and the space charge parameter C is given by 


C 


IT 
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(4) 
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As previously mentioned, equation (3) is solved by Picard iteration. The particular 
method, however, was that of Clenshaw and Norton (ref. 2), This method is appropriate 
here because it allows full use to be made of the many desirable features of Chebyshev 
expansions. 

An expansion in Chebyshev polynomials enables one to minimize the number of data 
points at which n(x) must be evaluated for a given average error. This is particularly 
important in the problem because n(x) is obtained from a Monte Carlo calculation, as 
described in the next section. 

In reference 1, the Chebyshev expansion n(x) = Saj^Tj^(x) was obtained and then 
transformed to a power series expansion in x before equation (3) was integrated. For 
any reasonably high degree polynomial, however, the coefficients of the power series 
expansion are badly conditioned; that is, they become so large that all precision is soon 
lost due to machine truncation errors. On the other hand, the coefficients aj^ of the 
Chebyshev series expansion are well behaved. It can be shown that [aj^l < 2M where M 
is the maximum value of n(x). Since the Chebyshev polynomials Tj^(x) are such that 
|Tk(x) I < 1, the error of truncating the series after a finite number of terms is readily 
apparent from the magnitudes of the first coefficients neglected. 

The power of the Clenshaw-Norton method is in its ability to employ only Chebyshev 
expansions and to eliminate the need for transforming to a power series expansion with 
the attendant loss of accuracy. Another desirable feature is that the solution (in this case 
cp(x)) is also given in terms of a Chebyshev expansion instead of the usual tabulation of 
numerical values. 


CONVERGENCE OF SOLUTION 

The convergence of interest here is that of the sequence (p^(x), . . . , <^’j^(x) 

derived from the Picard iteration of equation (3). If n(x) were a given function (not a 
stochastic function) , then, under the physical constraints imposed on the electron den- 
sity, a solution exists to the initial value problem. This knowledge gives no information, 
however, on the rate of convergence or the effect of stochastic variations of n(x) on the 
sequence { ^ effect of stochastic fluctuations on the convergence is discussed 

in appendix C. 

In practice, after about three iterations (depending, of course, on initial distribution 
<p^(x)) the potential distribution "settled down." Thereafter (succeeding iterations), the 
distribution fluctuated within a relatively small region of the function space defined by the 
solutions of equation (3). Successive iterations were then treated as independent trials. 
Their sample mean was accepted as the solution. 
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MONTE CARLO EVALUATION OF ELECTRON DENSITY AND CURRENT 


The general concepts of Monte Carlo calculations for electron-diode problems are 
discussed in reference 1. Basic to the procedure is to sample a great number of elec- 
trons emitted according to a given velocity distribution and obtain their averaged contri- 
bution to the electron density distribution and current. From these averages, an esti- 
mate of the desired diode characteristics is obtained. 


Initial Velocity Components 

At the beginning of each electron history, the initial velocity components are chosen 
in accord with a half-Maxwellian velocity distribution at the emitter. It is necessary, at 
this point, to mention that, in the sampling process, the histories of electron fluxes and 
not elements of electron density are being traced. The elements of electron density are 
not spatially invariant entities while the electron fluxes are. It is for this reason that the 
initial velocity components are not chosen from equation (2) but from the velocity distri- 
bution of electron flux, which in dimensionless cylindrical coordinates is 

2 2 

g(u,V)du dV = 4uVe"^^ ^ du dV (5) 

where u is the dimensionless velocity component parallel to the x-direction, and V is 
the transverse component. The marginal distribution functions of the random variables 
u and V are 


Gjj(t)dt = Gy(t)dt = 2te"^ dt 


( 6 ) 


Hence, the initial velocity components are chosen by the equations 



j 


where and Ry are, ideally, random numbers uniformly distributed over the range 
from 0 to 1 (see ref. 1, p. 23). Computer programs for choosing randomly from the 
range 0 to 1 result in, at best, sequences of pseudorandom numbers. These numbers, 
for the most part, however, are sufficiently random for Monte Carlo applications. For 
further information on pseudorandom numbers see references 3 and 4, 
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Distance to Collision 


Given an energy-independent mean free path X, the probability that an electron will 
travel a distance I without colliding is Hence, the distance to collision is 

obtained from the equation 


= -X In (8) 

where is again a random number between 0 and 1 (see ref. 1, p. 12). 

In this code, however, energy -dependent free path lengths X(E) are considered as 
follows: The interelectrode space is subdivided into a series of virtual cells separated 
by imaginary planes parallel to the electrodes. An average energy E is defined as the 
average kinetic energy the test electron would have in a cell assuming no collisions. On 
the basis of this average kinetic energy, a mean free path X = X(E) is obtained and em- 
ployed in equation (8) to obtain the distance to collision I in the cell. If I is greater 
than the distance along the electron trajectory to the cell boundary, the process is re- 
peated in the new cell (unless an electrode is reached). If is less than this distance, 
a simple search routine is employed to determine the correct corresponding to the 
point of collision. 


Angle of Scatter 

ENEC considers only electron-neutral elastic collisions; hence, the simplifying as- 
sumptions of infinite-mass (stationary) target particles is employed. Since the scatter- 
ing is anisotropic, however, the angle of scatter is a function of the angle of incidence. 
The polar angle of incidence 0^ = tan”^(V/u) is known at each collision. Because of the 
symmetry of the configuration, the azimuthal angle of incidence may be arbitrarily as- 
sumed to be = 0. 

It is, of course, assumed that the differential cross section (t( 0,E) is known or can 
be approximated. The methods used to smooth, fit, and tabulate the data for input to this 
code are discussed in the first section of ENEC CODE (also see appendix B). 

The cumulative distribution function of the random variable 0 can be defined, know- 
ing ct( 0,E), as 


where 


P[0 < 9] 



sin 0’ d0’ 

a(E) 


( 9 ) 


6 



r 


a 


/'IT 

(E) = J a(0,E)sin 9 d9 

0 


( 10 ) 


In the present case, it was decided to sample cos 9 instead of 9. From equation (9), 
the cumulative distribution function of cos 9 may be obtained by a simple transformation: 


P[cos 0 < t] 


■fr- 


(cos~n*, E) 


a(E) 


dt’ 


As further explained in the section Angular d istribution of scatter in ENEC CODE, the 
random variable cos 0 is tabulated as a function of a uniform distribution R between 
-1 and 1 by solving the equation 



a(cos~n,E) .. 
a(E) 


( 11 ) 


for cos 0 (see eq. (A16) of ref. 1 and the accompanying discussion). The azimuthal 
angle of scatter cp is chosen from a uniform distribution over the range from 0 to 2ir. 

With the scattering angles after collision about the incident direction denoted by 
primes, the transformation back to the coordinate system of figure 1 gives the following 
expression for the cosine of the polar scattering angle: 

cos 0 = cos 0Q cos 0' - sin 0^ cos cp' Vl - cos^0' (12) 


Sampling from Electron Histories 

Density . - The locations of the data points Xj^, where the electron density is sam- 
pled, are specified by the CHEBY subroutine. The contribution of a test electron (unit of 
flux) to the density at a given Xj^ is 


1 

TT^/^u(Xjj) 


(13) 


for each passage past Xj^, where 


(14) 


u(x) + cp{^) - 
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X is the position of the last event (collision or emission), and u is the initial velocity 
o 1/2 o 

of this trajectory at x^. The n ' results from the nondimensionalyzing factor em- 
ployed for the velocity (see definitions of u and V in appendix A); note that equa- 
tion (13), averaged over the initial distribution g(u,V) (eq. (5)), gives n(0) = 1, as 
assumed. 

The sample density at Xj^ at the end of one iteration (N^ histories) is given by 


n(x^) = 




"l<\> 


(15) 


where the sum over i (flux passages past Xj^) may be greater than, equal to, or less than 
Nq because of collisions and/or reflections from the potential field. 

After each iteration, the density distribution n(x) is obtained by fitting a Chebyshev 
ejqjansion to the sampled values of n(x^) . 

Current to coll^^tor. - The ratio of current density to the collector to the emitted 
current density J/Jq is computed at the end of each iteration from 



(16) 


where is the number of test electron fluxes reaching the collector. 

Collisions and flux passage^. - The total number of collisions is tallied in subrou- 
tines XIC and XICTP for each iteration, and the sample means are presented in the out- 
put (see sample problem, appendix F). 

The total number of flux passages at each data point is similarly tallied and pre- 
sented. These data have proved helpful from both a heuristic point of view and for de- 
bugging. 


General Programing Features 

While the program details are described in the section ENEC CODE, clarification of 
certain features may help to connect the analytical description with the code logically. 

Equations (7) and (8), which define the random variables u , v , and are not 
used in the program functions. Instead, the function -In R was tabulated for 1024 equi- 
distant values of R (0 < R < 1) . The application of a table look-up is five times as fast 
as evaluating -In R on the Lewis IBM 7094 n. While the coarse graining of the initial 
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velocity distribution is not deleterious in the present class of problems, care must be 
taken if this method is used to treat, for instance, inelastic effects, because of the trun- 
cation of the Maxwellian tail. For the same reason of improving computing efficiency, X 
was tabulated as a function of E, and cos 9 was tabulated as a function of R and E. 

Instead of standardizing on one formula, several different quadrature formulas have 
been employed for increased computational speed. The types of formulas employed were 
dictated^ in part, by the desirability of tabulating the potential distribution <p(x) at pre- 
determined sets of mesh points before each iteration, instead of evaluating the potential 
each time from its Chebyshev expansion. 


VELOCITY AND ENERGY DISTRIBUTION HISTOGRAMS 

After obtaining a solution for <p(x), the program was rerun for the express purpose 
of sampling the distribution functions. Rerunning did not involve any further iterations, 
but simply followed enough electron histories to obtain sufficient statistics. 

The histograms were obtained by first accumulating a sample of 250 histories of the 
random variable. This sample was then ordered, and from the resulting empirical dis- 
tribution the boundaries of 10 cells were chosen on the basis of equal probability. Then 
as many additional histories were tallied in the deciles as computer time and prudence 
would allow. In any event, this procedure precluded choosing cells wherein the subse- 
quent sample was too small. 

The original histograms obtained by sampling had to be modified. Whereas the veloc- 
ity distribution of the electron flux emitted at x = 0 is sampled, the histogram shows the 
velocity distribution of the electrons. The same is true for the energy distribution. 

More explicitly, a sample may be taken from a flux distribution function g(u,V,x); but 
this distribution function is related to the density distribution function f(u, V,x) by 

g(u, V, x) = Auf(u, V, x) (17) 

where A is a normalization factor. If, at given x, u is independent of the other veloc- 
ity components, marginal distributions of u become 

g(u,x) = Auf(u,x) (18) 


or 


f(u,x)o=Sfili2l 


u 


(19) 


Hence, to obtain the histogram of f(u,x), it is only necessary to divide the ordinate in 
each cell of g(u,x) by u, evaluated at the center of the cell, and to normalize. 

If, at given x, the kinetic energy is independent of polar and azimuthal angles. 


/ 2it /•?r/2 

dcp / uf(E, 0,<p,x)sin 0 d0 
u 0 


= A f duip f ^ cos 0 f(E, 0,^,x) sin 0 d0 

*'0 


= A’E^/^f(E,x) 


( 20 ) 


( 21 ) 


or 


f(E,x)oci®2^ (22) 

The procedure for obtaining the histogram of f(E,x) is then directly analogous to that 
for f(u,x). 
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ENEC CODE 


GENERAL FEATURES 
Cross Sections 

The user must supply either a functional expression or data for the differential 
scattering cross section a(0,E). In the case of data, a functional expression is first ob- 
tained by the method of spline interpolation, as described in appendix B. 


Geometry, Subdivisions, and Functional Tabulations 

ENEC is programed for a one -dimensional geometry as depicted in figure 2. Four 
sets of subdivisions are employed: 

Set a: 1024 Equally spaced subdivisions 
Set b: NS equally spaced cells 

Set c: Nonequidistant subdivisions defined by the abscissas of the three-point 
Gaussian quadrature formula (see G)uadrature Formulas) in each cel l 
Set d: Nonequidistant subdivisions defined by N1 arguments (called XD) of the 
Chebyshev curve fit (see Quadrature Formulas) 

The potential distribution ^(x) is tabulated at the 1025 equally spaced boundaries of 
set a; ^(x) is also tabulated on points defining sets c and d. These tabulated values of 
(p(x) are employed in the trajectory calculations (see Quadrature Formulas). 
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The cells, set b, are employed to obtain a spatial average of the electron kinetic 
energy (see section Distance to Collision) . Note that set c is really a set of subsets of 
set b (one subset per cell) . In addition, the potential minimum (PHIMIN) and its location 
(XMIN) are tabulated for each cell . The boundaries of set b are tabulated in the array 
XB. Sets c and d are also depicted in figure 2 for NS = 8 and N1 = 11. 


Trajectories 

In figure 3, the possible trajectories in a cell are shown for an electron moving to- 
ward the collector in a retarding potential field (fig, 3(a)). In figure 3(b), the electron is 
just passing into the cell and has an initial location XO, the cell boundary, while in fig- 
ure 3(c), the trajectories begin at the location of the last collision. In figures 3(b-3), 
3(b-5), and 3(c-3) to 3(c-6), the square of the u-component of velocity USQ is such that 
a turning point XTP occurs in the cell. 


I 
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I 


I 


□ Denotes collision 
O Denotes previous collision 




I i I 
XO XC XTP 

(b-5) 


I 

XO 


(c-1) 


(b) Trajectories originating at cell boundary. 

I 

f 



I 

XO 

(c-2) 


I 

XC 




(c-3) 


(c) Trajectories originating at collision in cell. 



I I I 
XO XC XTP 

(C-5) 



Figure 3. - Possible trajectories in a cell for electrons moving toward the collector. Potential, (p; initial cell boundary, XO; collision 
coordinate, XC; turning point coordinate, XTP; initial and final coordinates for numerical integrations, i and f, respectively. 
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Quadrature Formulas 


As mentioned in the ANALYSIS , several different quadrature formulas are employed 
to compute distances along the trajectories in order to reduce computing time. To com- 
pute distances along a trajectory, the integrand of equation (23) must be evaluated at 
every argument of the quadrature formula employed: 


1=1 y 1 + dx (23) 

J 1 u^ - (p{x^) + <p(x) 

Hence, computing time can be reduced by minimizing the number of arguments needed for 
a given accuracy and by having (p(x) tabulated at all such arguments. 

Gauss’s Quadrature Formula (ref. 5, p. 150) is one of the most efficient known for 
the numerical integration of a well-behaved function. The abscissas, however, are irra- 
tional numbers (in general). This fact together with the desireability of using only tabu- 
lated values of <p(x) precludes the extensive use of Gauss's formula here. This formula 
(QUAD) is used in ENEC, therefore, only for trajectories of the type depicted in figure 
3(b-l) and 3(b-2). 

For trajectories such as those depicted in figure 3(c-l), Simpson's Rule (ref. 5, 
p. 137) (QUADS) is employed between the limits of integration i and f over an equidis- 
tant subset of set a (see section Geometry, Subdivisions, and Functional Tabulation). 

This includes trajectories such as those depicted in figure 3(c-3) between the limits ig 
and f 2 . 

Simpson's Rule is not appropriate, however, for trajectories wherein a turning point 
exists (or would exist in the absence of collisions), for then the integrand (eq. (23)) be- 
comes infinite at the turning point (XTP) . It can be shown that, in the neighborhood of a 
turning point, the denominator in equation (23) goes to zero as (x - XTP)^''^, Hence, a 
Newton-Cotes type x^"^^ -weighted quadrature formula (QUADTP) was derived (see 
ref. 6), which did not require the integrand to be evaluated at x = 0. 


Tabulated Distributions 

Exponential distribution . - The exponential distribution is used to obtain the initial 
velocity components and the distance to collision, (eqs. (7) and (8)). Instead of evaluating 
-ALOG(R) for every random number R generated, -ALOG(R) is tabulated in array VEL 
for 1024 equidistant values (0 < R < 1) . 
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Free path length . - The energy-dependent mean free path X(E) is obtained from the 
total collision cross section a(E) (eq. (10)), by the formula X(E) = 1. 01/a(E). Then 
X(E) is tabulated in array MFP over the required range of E. 

Angular distribut ion of scatter. - The cosine of the scattering angle (cos 6) is ob- 
tained at a mesh of points (R, E) by numerical solution of equation (11) A surface fit of 
cos 9 = f(R, E) is then obtained (see appendix B), and subsequently cos 9 is tabulated in 
the two-dimensional array DISTB(R, E). 


Initial Velocities 

The square of the initial velocity components USQ (equal to USQO at x = 0) and VSQ 
is exponentially distributed (see eq. (7)). Use of the array VEL (see section Exponential 
distribution) is made in the following way: 

(1) Choose a random number R. 

(2) Let I = [1024.* R + 1. 5], where the brackets imply a truncation to the nearest 
integer. 

(3) Then USQO = VEL(I). 

(4) Repeat steps (1), (2), and (3) for VSQ. 


Elastic Collisions 

Distance to collision . - After determination of the average electron kinetic energy E 
in a given cell, or between the last collision and cell boundary (see section Distance to 
Collision), a distance to collision FPATH is chosen in the following way: 

(1) An integer I associated with E is obtained. (The value of this integer I de- 
pends on how X(E) is tabulated as a function of E in array MFP.) 

(2) Choose a random number R. 

(3) Let J = [1024.* R + 1. 5] . 

(4) Then FPATH = MFP (1)* VEL (J). 

Location of collision . - If the distance to collision (along the trajectory) is less than 
the distance to a cell boundary (along a trajectory in the absence of collisions) as, for ex- 
ample, depicted in figures 2(b-2) and 2(c-2), the collision location XC must be deter- 
mined. A simple binary search that uses either Simpson's Rule (XIC) or a x”^'^^- 
weighted Newton-Cotes formula (XICTP) is employed to obtain XC. In the neighborhood 
of a turning point, XICTP is used. 
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Angle of scatter . - After locating XC, the angles of scatter relative to the electron 
velocity vector before scatter are obtained in the following way: 

(1) A random number R is chosen. 

(2) The azimuthal angle (p' is given by 2-nB.. 

(3) The total kinetic energy E(XC) of the electron at XC is computed. 

(4) Another random number R is chosen. 

(5) Integers I and J are associated with R andE(XC), respectively. 

(6) Cos 0’ is given by DISTB(I, J). 

The cosine of the scattering (polar) angle (cos 0) in the original system of coordinates is 
then obtained from equation (12) where 


u(XC) 
cos 0 = ^ 

^E(XC) 

V 

sin 0Q = 

^E(XC) 

and u(XC) and V are the velocity components before collision. 


Generation of Electron Histories 

During the complete trajectory of the electron from the emitter to the collector, or 
back to the emitter, the following items are tallied: 

(1) The contribution to the electron density (see section Density) for each passage 

past an argument XD of the Chebyshev curve fit (see Solution of Differential 
Equation) 

(2) The number of electrons NTHRU reaching the collector 

(3) The total number of collisions in one iteration 


Random Number Generation and Selection 

The method employed at Lewis for generation of the pseudorandom number R on the 
unit interval is of the "congruential multiplicative" type (ref. 3). Most computing in- 
stallations have their own library routine for this purpose, but for completeness, and 
possibly for those users who would prefer not to change the calling sequence, the basic 
machine language (MAP) program RANDOM used in ENEC is given in appendix E . 
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Solution of Differential Equation 


After each iteration, the CHEBY routine is called to obtain a Chebyshev polynomial 
approximation (ref. 7) for the electron density distribution DEN(X). Then DEN(X) is em- 
ployed in the program CLN2S to obtain a new potential distribution by the method of 
Clenshaw-Norton (see CONVERGENCE OF SOLUTION) . The iterations are continued af- 
ter convergence (see VELOCITY AND ENERGY DISTRIBUTION fflSTOGRAMS) for pur- 
poses of averaging. 


Averaging Iterations 

As discussed in the section VELOCITY AND ENERGY DISTRIBUTION fflSTOGRAMS, 
convergence is only achieved in a stocastic sense dependent on the random error asso- 
ciated with the Monte Carlo evaluation of the density (see appendix C). Hence, in CLN2S, 
as soon as the variation in Chebyshev coefficients for <j?(x) falls within a precalculated 
range for two successive iterations, "convergence" is assumed. At this point KI more 
iterations are performed, the resulting densities at each XD are averaged, and a final 
evaluation is made to obtain (p{x). In addition, the KI sample collector potentials and 
currents are averaged and their standard deviations are computed. 


DIRECTIONS FOR ENEC USERS 
Preparation of Input Tables 

Three tables must be constructed and punched on cards in a binary format before 
ENEC may be used. These tables contain information needed by ENEC to compute the 
initial velocity components of the electrons, the distance to a collision (free path length), 
and the angle of scatter after a collision (see Tabulated Distributions). Listings and de- 
scriptions of the programs that construct these tables (CVEL, MFP, ARGON, and 
ARGINV) are given in the section Auxiliary FORTRAN IV Program Descriptions. The 
SPLINE curve and surface fitting subroutines (SPLINE and SPLIN2) needed for interpola- 
tion, and the machine language subroutine for punching binary formatted cards (BCDUMP) 
are given in appendixes B and E, respectively. Samples of input and output for programs 
CVEL, MFP, ARGON, and ARGINV are given in appendix F. 

Exponential distribution. - To obtain the table of the exponential distribution (see 
section on Exponential Distribution) on punched cards, load and execute the programs 
shown in figure 4. The deck of output cards is to be loaded with ENEC, as shown in 
figure 5. 
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Figure 4. - Deck configuration for 
exponential distribution table. 


Free path lengths . - To obtain the table of X(E) (see section Free Path Length, p. 14), 

load and execute the programs shown in figure 6. The input consists of known values of 
1/2 

CT (E ' ) for the gas under consideration. The deck of output cards is to be loaded with 
ENEC. 

Angular distribution of scatter . - Two steps are involved in obtaining the table of 
the angular distribution noted in the section Angular Distribution of Scatter. The first 
step is to execute the deck configuration of figure 7 (a) . The input to this deck is a table 
of known values for the differential cross section for the gas under consideration ct( 0,E) 
(see section Angle of Scatter) . The output from this deck (values on the surface defined 
by eq. (11)) is then input to the deck configuration shown in figure 7(b). This second deck 
configuration produces the table of angular distribution on punched cards that is to be 
loaded with ENEC. 


Input Data - Problem Preparation 

Input to ENEC consists of two parts. In the first part, three tables are constructed 
and punched on cards by the deck configurations of the preceding section. The binary 
formatted cards punched by BCDUMP are read by subroutine BCREAD given in appen- 
dix E. 

The second part of the input to ENEC consists of variables that define the particular 
diode configuration to be studied. This input is formatted by the FORTRAN NAMELIST 
statement. 
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Figure 6. - Deck configuration for table 
of energy-dependent mean free path. 



(a) First step. (b) Second step. 


Figure 7. - Deck configuration for obtaining angular distribution table. 
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The following outline gives the order (see fig. 5), format, and description of the in- 
put data of ENEC: 

(1) BCREAD (VEL(1),VEL(1024)) 

The table of values computed by program CVEL 

(2) BCREAD (MFP(1),MFP(126)) 

The table of values computed by program MFP 

(3) BCREAD (DISTB(1, 1), DISTB(64, 64)) 

The two-dimensional table of values computed by programs ARGON and ARGINV 

(4) NAMEUST/INI/NO, N1,KI, ALPHA, CONST, NFLAG, BC, KODE, MODE, ERROR, 

NS, KFLAG, TEMPK, LFLAG 

NO Number of particles to be processed for each iteration 

N1 Number of Chebyshev data points to be used; N1 < 17 and odd 

KI Number of iterations after convergence in CLN2S to gather statistics 

on various parameters; KI < 20 

ALPHA Ratio of mean free path to electrode spacing when KFLAG = 1; P^L 
when KFLAG = 2 

CONST Constant C in Poisson’s equation 

NFLAG Control on reading of AI array, NFLAG = 1; AI is read by 
NAMELIST/IN2/. NFLAG = 2; AI is read by BCREAD. 

BC Boundary condition on differential equation, ^’(0) = BC when 

MODE = 1, and <p(l) = BC when MODE = 2 

KODE Control on output from CLN2S (see section ENEC output). KODE = 1; 

CLN2S plots ^(x) and writes Chebyshev coefficients for each itera- 
tion during convergence. KODE = 2; no output from CLN2S. 

KODE = 3; no plots, but Chebyshev coefficients are written for each 
iteration during convergence. 

MODE Control on boundary condition BC. MODE = 1, cp'(0) = BC; 

MODE = 2, (^(1) = BC. 

ERROR Convergence is assumed in CLN2S when the maximum difference be- 
tween the corresponding coefficients AI for two successive itera- 
tions is less than ERROR. 

NS Number of cells; must be 1, 2, 4, 8, or 16 

KFLAG Control on the use of the MFP table. KFLAG = 1; constant ALPHA 

assumed. KGLAG = 2; ALPHA dependent on MFP table. 

TEMPK Temperature constant 

LFLAG Control on the calculation of the scattering angle; LFLAG = 1, iso- 
tropic scattering angle assumed; LFLAG = 2, scattering distribu- 
tion used, DISTB. 
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(5) The initial coefficients for the Chebyshev fit of ^(x) are read at this point. 

NFLAG controls the format of this input. 

For NFLAG = 1, the input format is: 

NAMEIIST/IN2/AI 

For this case the data read into AI are usually the coefficients of a 
Chebyshev fit for a straight line satisfying the boundary condition BC. 

For NFLAG = 2, the input format is: 

FORMAT (12A6) 

BCREAD (AI(1),AI(20)) 

BCREAD (DA(1),DA(20)) 

BCREAD (B(1),B(20)) 

The first card of this input contains 72 Hollerith characters describing the 
ENEC run producing the AI data. 

The next three binary cards contain the coefficients of the Chebyshev fit of 
<p(x), (p'{x), and (p'*(x) produced by a previous ENEC run. Although the coeffi- 
cients for (p'(x) and cp"(x) are not needed, they are read in to keep all coeffi- 
cients generated by a particular ENEC run together. 

(6) For multiple runs, the NAMELIST statement /INI/ is read until the input is ex- 

hausted. 


ENEC Deck Configuration 

The deck and input setup for execution of ENEC is shown in figure 5 (p. 18). Deck 
RANDOM is the random number generator given in appendix E, and deck FSQR is a fast 
square root subprogram designed specifically for Monte Carlo work, A discussion of 
square root subroutines is given in appendix D, while the square root subroutine itself is 
presented in appendix E. The printer plot subroutines are used to give plots of the func- 
tions <p(x), <p'(x), and <p”(x) computed by ENEC. Reference 8 presents the necessary 
subroutines. The user may, if he wishes, remove all calls to PLOTF from ENEC thus 
eliminating the need for the plot subroutines. The input has been described in the pre- 
vious section and must be placed in the order shown in figure 5. 

ENEC Output 

The output generated by ENEC consists of printed listings and punched cards. The 
printed output lists the input read by NAMEUST/INl/ and the initial coefficients of the 
Chebyshev fit of <^(x) along with the following information computed by the code: 
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(1) The N1 Chebyshev data points are listed with the mean and standard deviations 
of the number of particles crossing each data point. 

(2) The mean and standard deviations of the N1 coefficients of the Chebyshev fit of 
^(x), <p’(x), and ^”(x) are listed. 

(3) The mean and standard deviations of <^(x), <p’(x), and cp"{x) are listed at the 
Chebyshev data points. 

(4) The mean and standard deviations of the current, the voltage [<??(!)] , <p'(0), and 
KNTR (where KNTR is the number of collisions of the particles per iteration) are listed. 

(5) Values of . and x . for the mean (p(x) are listed. 

mm tnxn 

(6) Values of <p(x), cp'ix), and n(x) [<p”(x)/C] are listed at 21 equally spaced points 
over the range from 0 to 1. 

(7) Printer plots (ref. 8) are given for the three functions (p(x), cp'ix), and <p"(x) 
over the range from 0 to 1. 

Depending on the input variable KODE, output may be obtained from subroutine CLN2S. 
For each iteration during convergence in CLN2S, the user may have printed either the 
coefficients of the Chebyshev fit of (p(x), cp'ix), and cp"ix) and printer plots of these 
functions, or just the coefficients of these functions. Thus, the user may ’’see" the con- 
vergence process take place in CLN2S. 

The punched output consists of three binary cards, punched by BCDUMP, containing 
the mean coefficients of the Chebyshev fits of cpix), cp'ix), and cp"ix). The information 
contained on these cards is useful for input to other ENEC runs and for further investiga- 
tion into the properties of a particular diode by other computer programs. For an exam- 
ple of the printed output, see appendix F. 



0 4 8 12 16 20 24 28 32 36 


Number of collisions experienced by an electron 
Figures. - Execution time. 
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Execution Time 


An estimate of the execution time for ENEC is difficult to predict because the number 
of parameters used to define the diode and Monte Carlo solution is large. Figure 8 does 
show, however, that the execution time for one iteration does vary linearly with respect 
to the number of collisions that the electrons encounter. Experience with ENEC has 
shown that an average run usually takes about 13 iterations: 3 for convergence and 10 to 
gather statistics. The number of collisions an electron will make is, again, difficult to 
determine. Generally, the number of collisions will increase when the voltage is in- 
creased or when the mean free path is decreased. The digital computer used to run ENEC 
was an IBM 7094 11-7044 direct couple system. Therefore, the execution times are based 
on this system. 


PROGRAM DETAILS 
ENEC Labeled COMMON 

A description of all FORTRAN variables appearing in labeled COMMON in ENEC is 
given in table I. The COMMON label is listed with the variables in the order and with the 
dimension information used in ENEC. The cross reference between ENEC programs and 
labeled COMMON are shown in table II . 
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TABLE I. - DESCraPTION OF FORTRAN VARIABLES APPEARING IN ENEC COMMON 


COMMON 

label 

FORTRAN 

variable 

Description 

/CMAIN/ 

NO 

Number of particles to be processed for one iteration 


NI 

Number of Chebyshev data points 


ALPHA 

Ratio of mean free path to electrode spacing or the 
product P^L 


CONST 

Constant in Poisson^s equation 


NS 

Number of cells 


AI(20) 

Array of coefficients of Chebyshev fit to <p(x) 


VEL(1024) 

Array of values of -In R, where R is equally spaced 
between 0 and 1 (see section Exponential distribution) 


NK 

Counter of iterations after convergence 


MFP(126) 

Array of values of X(E) (see section Exponential 
distribution) 


KFLAG 

Control used in choosing mean free path; KFLAG = 1 
for constant ALPHA; KFLAG = 2 for ALPHA de- 
pendent on MFP array 


THETA 

TEMPK/11600. 0 


DISTB(64, 64) 

Gas scattering angle distribution (see section Angular 
distribution of scatter) 


LFLAG 

Control used in choosing scattering angle; LFLAG = 1 
for an isotropic scattering angle; LFLAG = 2 for 
scattering angle chosen from DISTB 

/CITER/ 

N(20) 

Counter for the number of particles passing each 
Chebyshev data point 


Y(20) 

The tally count at each Chebyshev data point; this 
array contains the new cp^^(x) after ITER and de- 
termines the new Chebyshev fit 


KNTR 

Counter for the number of collisions for each iteration 


VSQ 

y2 


FPATH 

Mean free path 


IBO 

Index of cell boundary behind particle 


IBF 

Index of cell boundary ahead of particle 


SIGMA 

Control indicating direction of particle; SIGMA = 1 if 
the particle is traveling to the right; SIGMA = -1 if 
the particle is traveling to the left 
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TABLE I. - Concluded. DESCRIPTION OF FORTRAN VARIABLES APPEARING 

IN ENEC COMMON 


COMMON 

label 

FORTRAN 

variable 

Description 

/CITER/ 

ID 

Index of next Chebyshev data point ahead of particle 


CRRNT 

Ratio of number of particles processed to number es- 
caping to the right of the diode per iteration 

/CMNPHI/ 

NDEL 

1024/NS; number of equally spaced data points in each 
cell 


DELX 

1./1024; distance between equally spaced data points 


XB(20) 

Cell boundaries 


TPHIQ{3, 20) 

Functional values of (p(x) at three Gaussian data points 
for each cell 


IMIN(20) 

Index of location of for each cell 

^min 


PfflMIN(20) 

for oach cell 

^min 


XMIN(20) 

Location of for each cell 


TPHID(20) 

Functional values of (p (x) at the Chebyshev data points 


TPmX(1026) 

Functional values of (^(x) at the equally spaced data 
points 

/CCLN2S/ 

DA (20) 

Array of coefficients of Chebyshev fit to <p^(x) 

/CCHEBY/ 

B(20) 

Array of coefficients of Chebyshev fit to 


XD(20) 

Chebyshev data points 

/CCELL/ 

NCELL 

Cell number particle is in 


XO 

Position of particle 


XF 

Position of cell boundary ahead of particle 


lO 

Index of particle position 


IF 

Index of position of cell boundary ahead of particle 


USQO 

2 

'^o 

/CXITP/ 

ITP 

Index of position of turning point 


XTP 

Position of turning point 

/CXICTP/ 

IC 

Index of position of collision 


xc 

Position of collision 

/CPATH/ 

EV 

Energy of particle 
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TABLE n. - LABELED COMMON AND ENEC PROGRAM CROSS REFERENCE 


ENEC 

COMMON label 

program 

/CMAIN/ 


/CMNPHI/ 

/CCLN2S/ 

/CCHEBY/ 

/CCELL/ 

/CXITP/ 

/CXICTP/ 

/CPATH/ 

MAIN 

X 

X 

X 

X 

X 





ITER 

X 

X 

X 



X 




MINPffl 

X 


X 


X 





CELL 


X 

X 


X 

X 

X 

X 


STOSS 

X 

X 

X 



X 


X 

X 

PATH 

X 

X 

X 



X 



X 

XIC 


X 

X 





X 


XICTP 


X 

X 




X 



XITP 



X 



X 

X 



QUAD 

X 

X 

X 



X 




QUADTP 


X 

X 



X 




QUADS 

CLN2S 

X 

X 

X 

X 

X 

X 




CHEBY 

X 

X 


1 

X 





PHI 

X 









DPHI 

X 



X 






DENS 

DISCRl 

DISCR2 

PLOTF 

PLOTYX 

SORTYX 

SCALEY 

X 
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ENEC FORTRAN IV Program Descriptions, Flow Charts, and Listings 

Table HI presents a brief description of all FORTRAN programs used in the ENEC 
code. The MAP coded programs (FSQR, RANDOM, BCREAD, and BCDUMP) are pre- 
sented in appendix E. Printer plot programs used by PLOTYX are given in reference 8. 

A directed graph, of the programs given in table III is shown in figure 9. This figure de- 
picts the overall interrelation of the ENEC FORTRAN programs. Listings and flow 
charts (figs. 10 to 36) for the FORTRAN programs are presented in the following sections. 
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TABLE m. - DESCraPTION OF ENEC FORTRAN IV PROGRAMS 


Program 

Description 

MAIN PROGRAM 

Controls the flow of the Monte Carlo solution of planar electron- 
diode problems; reads all input and writes most output 

CLN2S 

Solves second-order ordinary differential equations of the form 
y*^(x) = f(x,y,y') 0 < x < 1 by the method of Clenshaw -Norton 
using Picard iteration 

ITER 

Initializes variables for the start of an iteration; chooses initial 
conditions of an electron entering the diode; determines the 
cell and cell boundaries that the particle is entering; counts 
the number of particles passing through the diode; normalizes 
the functional values of (p”{x) at the Chebyshev data points 

MINPHI 

Computes the cell boundaries and the functional values of cp{x) 
at the necessary data points; locates along with the 

corresponding data point and index for each cell 

CELL 

Follows the path of a particle through a cell 

STOSS 

Determines conditions of a particle after collision such as direc- 
tion and velocity 

PATH 

Computes the energy and mean free path of a particle 

XIC 

Determines the position of a collision when a turning point is 
not involved 

XICTP 

Determines the position of a collision when a turning point is 
involved 

XITP 

Determines the position of a turning point 

QUAD 

Computes the distance of path across a cell by applying a three - 
point Gaussian quadrature 

Q0ADTP 

Computes the distance of path between two data points where one 
of the data points is a turning point; a three -point open-ended 
quadrature is used near the turning point, while Simpson's 
Rule is used otherwise 

QUADS 

Computes the distance of path between two data points by apply- 
ing Simpson's Rule; neither data point is a turning point 

CHEBY 

Initially computes the Chebyshev data points for a given N1 and 
subsequently obtains coefficients for the Chebyshev fit of 
<p"(x) 

PHI 

Computes the functional values of cp (x) from the Chebyshev fit 

DPHI 

Computes the functional values of <p^(x) from the Chebyshev fit 

DENS 

Computes the functional values of (^”(x) from the Chebyshev fit 

DISCRl 

Averages and computes the standard deviation of values stored 
in a two-dimensional array 

DISCR2 

Averages and computes the standard deviation of values stored 
in a one-dimensional array 

PIX)TF 

Plots functions described in the calling vector on a single page 

PLOTYX 

Plots values stored in arrays on a single page 

SORTYX 

Sorts arrays to be plotted 

SCALEY 

Scales the ordinate array so that the plot will be contained on a 
single page 




Figure 9, Directed graph of FORTRAN programs used in ENEC code. 







































$I0FTC MAIN 


COMMON /CMAIN/ NO , N1 , ALPHA » CONST , NS , A I (20) ,VEL( 102^ ) , NK , MFP ( 126) 
« KFLAGf THETA, OIST0 (64,64) tLFLAG 
COMMON /CITER/N(20)f Y( 20) , KNTR, VSQ , FP ATH , I BO , I 8F , S IGMA , ID,CRRNT 
COMMON /CMNPHI/NOELfOELXfXB! 20) , TPH I Q ( 3 , 20 ) , I MI N( 20 ) , PH IMI N C 20 ) , 
«XMIN( 20), TPHID<20) ,TPHIX(1026) 

COMMON /CCLN2S/ 0A(20) 

COMMON /CCHFRY/ B(20),XD(20) 

DIMENSION MAK 20,2 0) ,MOA( 2 0,20) , MB (20,20) , STD A I (20) ,STD0A(20I , 

* STDR( 20) , MY (20, 20), MO Y( 20,20 ), MOD Y( 20,20) ,STDY(20) ,STD0Y(20) , 

* STOOD Y( 20) ,DY( 20) ,00Y( 20) ,CURRNT<20) ,H(12),HH(l2).Xl(2l),Yl(2l) 

* DYK 2l),00Y 1( 21) ,FN( 20,2 0) ,FKNTR(20) ,FN1 (20) ,ST0N(20) 

DATA H/4*6H , 6HA I DAT,6HA READ,6H BY NA , 6 HMEL I ST , 4=^6H 

DATA HH/4*6H ,6HAI DAT,6HA FR0M,6H PRECE , 6HEDI NG ,6HRUN , 

<'3«6H / 

DATA M/1025/ 

REAL MFANC,MEANV,MAl ,MDA ,MB,MDPHlOt MF P , MY , MOY , MOOY 
EXTERNAL P HI , OPHI , OENS 

NAMELIST / IN1/N0,N1, KI , AL PHA , CONST , NF LAG , BC , KODE , MODE , ERROR , NS , 
^J'KFLAG, TEMPK,LFLAG 
NAMELIST /IN2/ AI 
CALL SAND(RO) 

C READ IN INITIAL DATA 

CALL BCREAD(VEL( 1),VFLC 1024) ) 

CALL BCREADl MFP( 1 ) ,MFP( 126) ) 

CALL BCREAO( DISTBC 1 , 1 ) , 0 I S TB ( 64 , 64) ) 

READ ( 5, INI) 

WRITE (6,202) 

WRITE (6, INI) 

IF(NFLAG.EQ. 1) GO TO I 

READ (5,101) H 

CALL BCREADIAK l),AI(20)) 

CALL BCREAO(OA( l),OA( 20) ) 

CALL 8CREA0( B( I ) , 0( 20) ) 

GO TD 2 

1 READ ( 5, IN2) 

2 WRITE (6, lOl) H 

WRITE (6,201) (AM I ) ,1=1,20) 

NK = 0 

THETA = TEMPK/U600.0 
5 CALL CLNZStKODE, MODE, 3C, ERROR) 

NK = NK 1 
DO 6 1=1, N1 
FN(NK, I) = N( I ) 

MY(NK,I) = PHUXO( m 
MOY(NK, I ) = DPHi { XD( I ) ) 

MOOY( NK, I) = OENSI XD ( I ) ) 

MAUNK, I ) = AI ( I ) 

MDA(NK, I ) = DA ( I ) 



6 HBiNKt I) = B( I ) 

FKNTRCNK) = KNTR 
CURRNT(NK) = CRRNT 
IF(NK.LT.KII GO TO 5 

CALL DISCRl(NK,Nl,MY,Y,STOY) 

CALL OISCRl(NK,Nl,MOY,DY,STOOY» 

CALL OISCRl(NKfNl.MOOY,OOY,STDOOYI 
CALL DISCRKNK,Nl,MAI ,AI,STOAI ) 

CALL OISCRUNK,Nl,MDA,OA,STOOA) 

CALL 01SCR1(NK,N1,MB>B*ST0B) 

CALL DISCRl(NK»Nl,FNfFNl,STON> 

CALL OISCR2(NK,FKNTRfFKNTRl,STOKTR» 

CALL DISCR2(NK,CURRNT,MEANC,STDC» 

MEANV = Y(N1I 
STDV = STOYCNU 
MDPHIO = OY( l» 

STOPHI = STOOY(l) 

DO 7 1 = 1,21 

Xim = FLOAT! I-n/20.0 
YHI) = PHI(Xim) 

DYim = DPHKXKin 

7 DOYl(I) = DENSfXlU n/CDNST 
PHIM = 0.0 

XM = 0.0 

00 9 1=1, M 

X = 0ELX*FL0AT(I-1) 

PHII = PHKXI 
IFCPHII.GT.PHIMI GO TO 9 
PHIM = PHI I 
XM = X 
9 CONTINUE 

WRITE (6,205) ( I , XOU ) ,FN1 ( I) , STONI I ) , 1 = 1 , N1 ) 

WRITE (6,203) ( I , AI (I ) , STOA I( I) ,0A( I ) , STDOA( I ) , BU ) ,ST 0B( I) , 

♦1=1, Nl) 

WRITE 16,207) ( I, Y( I ) , STOY( I ) ,OY<I ) , STDOYt I ) , OOYC I) ,STODOY ( I ) , I = I, 
♦Nl) 

WRITE (6,204) MEANC, STOC, MEANV, STOV, MDPHIO, STOPHI ,FKNTR1,ST0KTR 
WRITE (6,209) PHIM,XM 

WRITE (6,208) ( X I ( I ) , Y I ( I ) ,0Y1 ( I ) , OOYl ( I ) , 1=1 , 21 ) 

CALL BCDUMP(AI( l),AI(20n 
CALL BC0UMP(0A(1),0A(20)) 

CALL 8C0UMP(BI 1),B(20) ) 

CALL PL0TF(21,0.0,1.0,PHl ) 

CALL PL0TF(2l,0.0,1.0,0PHI ) 

CALL PL0TF(21, 0.0, 1.0, DENS) 

READ (5, INI) 

WRITE (6,202) 

WRITE (6, INI) 

00 8 1=1,12 

8 H( I ) = HH( I) 

GO TO 2 

101 FORMAT ( 12A6) 

201 FORMAT ( 4HKA 1= , /, ( 6E 20. 8) ) 

202 FORMAT ( IHl) 

203 FORMAT ( IHK, 13X, 7HMEAN AI , 13X, 7HSTD. AI ,13X,7HMEAN DA , 13X,7HSTD. 0 
♦A, 13X,6HMEAN B,14X,6HST0. B , / , IHK, / , 1 IH ,I2,6F20.8)) 
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204 FORMAT ( IHK, 30X, 4HMEAN t 16X ,4HSTD. , / /, 1 3X , 7HCURRENT , 2F20. 8 , // , 13X , 4 
*HV0LT,3X,2F20.8,//t 13X,5HDPHI0,2X,2F20.8,//,13X,4HKNTR, 3X,2F20,3) 

205 FORMAT ( IHl, 16X, 2HX0 ,20X,6HMEAN N,l 9X, 6HST0. N,///,(IH ,I2,F20.8,2 
♦F20.3) ) 

206 FORMAT J1HK,2H K, lOX , IHY , I 8X.2H0Y , I 7X , 3H0D Y ,/ , ( I H ,I2»3F20.8)) 

207 FORMAT ( IHK , 14X, 6HMEAN Y,14X,6HST0. Y t 13X , 7HME AN OY » 1 3X , 7HST0. DY, 
♦ 12X,8HMEAN DDY , 12X, 8HST0. DO Y, / , IHK, / , 1 1 H ,I2,6F20.8)) 

208 FORMAT (20H1EQUALLY SPACED DATA , /// ,1 OX , IH X , 1 9 X ,6HPHI ( X ) , 1 4X , 7HDPH 
*n X),13X,4HN(X) ,///,(4F20.8» ) 

209 FORMAT { IHK, 12X, 7HPHI M [N= ,F 15. 8 , // , 1 3 X ,5HXMIN= , F 17 . 8 > 

END 
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Figure 11. - Flow chart for 






I 



subroutine CLN2S. 








$I8FTC CLN2S 


SUBROUTINE CLN2S ( KOOE . MODE , BC » ERROR ) 

COMMON /CCLN2S/ 0A120) 

COMMON /CCHE8Y/ 8I20),X0(20I 

COMMON /CMAIN/ NO, N , ALPHA , CONST ,NS , A I I 20) , VEL ( 1024 » ,NKi MF P( 126 ) , 
♦ KFLAG, THETA, 0IST8(64, 64) ,LFLAG 
DIMENSION A(20) 

DATA NIT/10/ 

EXTERNAL PHI 
IF(NK.NE.O) GO TO 1 
CALL CHEBYd) 

HBC=.5*BC 
I DO 99 IT=1,NIT 
CALL. ITER 
CALL CHE6Y(2) 

IFIKODE.EQ.l.ANO.NK.EQ.O) CALL PLOTF ( 21 , 0 .0, I . 0 ,PHI ) 

8(N+1)=0 
DA(N+1 )=0 
DO 40 K=2,N 
R=K-1 

40 DA(K) =(8(K-1)-81K+1) )/I4.*R) 

DO 41 K=3,N 

R=K-1 

41 A(K) = (0AIK-l)-0A(K + in/(4.*R) 

GO TO (50, 60), MODE 

50 U=0 
V=0 

DO 51 1=2, N 

51 V = V + DAI I )♦{-!. )**I 
DAI 1)= 2.*IBC«^V) 

Al 2) = ( DAI 1 )-OAI 3) )/4. 

DO 52 1=2, N 

52 U = U + Al I )*(-!. )**I 
A( 1)= 2.*U 

GO TO 70 

60 U=0 
V = 0 

DO 61 1 = 3, N, 2 

61 U=U ♦ A( I ) 

A( 1 )=2.4( HBC-U) 

DO 62 1=4, N, 2 

62 V=V + A( I ) 

AI2) = HBC - V 

DAI 1)=DA( 3) *■ 4.*A12) 

70 ERRMAX = 0 
DO 80 K=1,N 
ERR=ABSIAI(K)-AIK)) 

80 IF( ERR. GT. ERRMAX) ERRMAX=ERR 
IF (IT.EQ.NIT) GO TO 81 
IF(KOOE.EQ.2.0R.NK.NE.O) GO TO 87 

81 WRITE (6,82) 
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82 FnRMATllHL,2H K , 9X » 2H A I , 19X 1 1 HB , 19X f IH A, 19 X , 2HDA/ IHJ ) 
WRITE 16,84) { If A I (I ) ,B{I ) ,Am ,DA( I ) ,1 = 1, N) 

84 FORMATdH ,I2,4F2Q.8) 

WRITE <6,85) ERRMAX,IT 

85 FORMAT ( 8HKERRMAX = , E 1 8 . 0, 5X ,3H I T=, I 4 ) 

IF ( IT.EQ.NIT) RETURN 

87 IF(NK.NE.O) GO TO 98 

AMM = ABS(A(N)) + A8S<A(N-D) 

IF(AMM ♦GT. ERROR) GO TO 92 
GO TO 98 

92 ADI = ABS( AI ( I )-A( in + AB S < A I ( 2 ) -A ( 2 ) ) 

AD2 = ABS(AI(N-2)) + ABStAKN-Bn 
IF( A0l.GT.A02) GO TO 98 
IF(N.GE.17) GO TO 98 
N=N^-2 

CALL CHEBY(l) 

98 DO 100 K=1,N 
100 AI(K) = A(K) 

IF( NK .NE.O) RETURN 

IF( ERRM AX. LT. ERROR) RETURN 

99 CONTINUE 
RETURN 
END 




















$I8FTC ITER 


SUBROUTINE ITER 

COMMON /CITER/N(20)f Y(20) , KNTR , VSQ , E P A TH , I RO » IBF,SIGMAt ID»CRRNT 
common /CMAIN/ N0,N1 , alpha ,CONST,NS t at ( 20) ,VEL( 1024 ) , NK , MF P ( 126 )t 
« KFLAG, THETA, 0IST3C64,64) tLFLAG 
COMMON/CMNPHI /NOEL , DEL X,XB ( 20) , TPHIQ( 3,20) , I M I N ( 20 ), PH I MI N ( 20 ) , 

* XMIN( 20 ) , TPHI0( 20) ,TPHIX ( 1026) 

COMMON /CCELL/NCELL,XO,XF ,10, IF,USQn 
DATA SQRTP I/U 77245385/ 

NTHRU = 0 

on 10 1=1,20 

N( I )=0 

10 Y{ I ) = 0.0 
CALL MINPHI 
KNTR=0 

00 33 II=l,N0 
CALL RAND(R) 

J = IFIX( 1024. ’i'R ) + 1 
USQO = VEL ( J ) 

CALL KANO(K) 

J = IFIX( 1024. «R) + 1 
VSO = VEL(J) 

NCELL = 1 
IBO = 1 

IBF = 2 
SIGMA = 1.0 
10 = 1 

20 10 = NOEL<^( IBO-n + l 
IF = NOEL=^=( I8F-D + 1 
CALL CELL 

IF( IBF.EQ. 1) GO TO 33 
IF( IBF.EQ.NS+1) GO TO 32 
NCELL = NCELL ^ IFIX(SIGMA) 

IBCJ = IBF 

IBF = IBO + IF IX( SIGMA) 

GO TO 20 

32 NTHRU = NTHRU + 1 

33 CONTINUE 

DO 34 I=l,Nl 

34Y(I) = (Y( I>/IFLGAT(NO)=«^SQRTPI ) )«=CONST 
CRRNT = FLOAT( NTHRUI/FLOAT (NO) 

RETURN 

END 










SIBFTC MINPHI 


SUBROUTINE M INPHI 

COMMON /CM N PH I /NOEL f OF LX, XB ( 20) ,TPHIQ(3,?0) , I M I N ( 20 ) , PH I MI N ( 20 ) , 

« XM IN 1 20) , TPH I 0( 20) , TPHIXI 1026) 

COMMON /CMAIN/ (NO , Nl , ALPHA , CONST , NS , AI ( 20 ) ,VE L ( 1024 ), NK,MFP ( 126 ) , 
KFLAG, THETA, 01 STB (64,64) ,LFLAG 
COMMON /CCHEBY/ B(20),X0(20) 
data N / 1024/, FN/ 1024.0/ 
data AA/7.74596692E-1/ 

OELX = l.0/FL0AT(NS) 

NOEL = N/NS 
DO 3 1=1, NS 

3 XBCn = DELX*FLOAT( I-l ) 

XB(NS+l) = 1.0 
DO 5 J=1,NS 

XX ^ ( XB (J +1 )-XB ( J ) )«0. 5 
TPHIQ(1,J) = PHI ( XBIJ ) 4-XX« ( l.O'-AA) ) 

TPHIQ(2,J) = PHK (XB( J)+XB( J + 1) )«0.5) 

5 TPHIQ( 3,J) = PHI(XB(J)^-XX*C1.0+AA) ) 

DO 1 1=1, Nl 

1 TPHIDC I ) = PHI (X0( I ) ) 

PHIMIN(l) = 1•OE^■30 
OELX = 1.0/FN 

M = N + 1 
J = 1 

DO 2 1=1, M 
XX=0FLX={^FL0AT( I- 1 ) 

IF(XX.LE.XB( J+1) ) GO TO 11 
J = J + 1 

PHIMIN(J ) = TPHIX( I-l) 

IMIN(J) = I“1 

XMIN(J) = OELX^FLOATC I~2) 

1 I TPHIX ( 1 ) = PHI ( XX ) 

IF(PHIMIN( J ) .LT.TPHIXC I) ) GO TO 2 
PHIMIN( J ) = TPHIX( I ) 

IM IN( J ) = I 
XMIN(J) = XX 

2 CONTINUE 
RETURN 
END 



XO = XBtIBO) 
XF - XB(IBF) 



Figure 14. - Flow chart for 
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subroutine CELL. 












$£BFTC CELL 


SUBROUTINE CELL 

COMiMDN /CCELL/ NC ELL t XO, XF , I U , I F t USQO 

COMMON /ClTER/N(20),Y(;?0),KNTR,VSQ,FPATHf IBO,IBF,SIGMA,ID,CRRNT 
COMMON /CMNPHI /NOEL »DE LX, XB ( 20 ), TPH I Q { 3 ,20 ), I MI N ( 20 ), PH I MI N ( 20 ) , 
^5^ XMIN(20) ,TPHI0(20),TPHIX( 1026) 

COMMON /cxrip/ ITP,XTP 
COMMON /CXICTP/ IC,XC 
COMMON /CCHERY/ B(20),X0(20) 

C ENTER CELL 

XO = X6i IBU) 

XF = X0( IBF) 

IFIUSQO •LT.-PHIMIN(NCELL) ) GO TO 71 
CALL PATH( 10, IF) 

3 S = QUAD( IZllll) 

30 IF(FPATH.LT,S) GO TO 4 

31 T = SIGMA^{XF-XO( 10) } 

IF( T.LT*0.0) GO TO 32 

TALLY = l*0/SQKT( USOO + TPHIOnO) ) 

Y( 10) = Y{ ID) + TALLY 
N( ID) = N( 10 ) + 1 
10=10+ IF IX( SIGMA ) 

GO TO 31 

32 FPATH = FPATH - S 
RETURN 

4 CALL X IC( 10, IF ) 

41 T = SIGMA*IXC-XO( 10) ) 

IF(T.LT*0.0) GO TO 5 

TALLY = 1.0/SQRT(USQ0+TPHI0(I0) ) 

Y( ID) = Y( ID ) + TALLY 
N{ ID) = N( ID) + I 
10=10+ IF IX( SIGMA) 

GO TO 41 

5 CALL STOSS 

TAU = XMIN(NCELL )-XC 

6 IF( TAU^SIGMA.LT.0.0) GO TO 9 

IF( USQO*LT.-PHIMIN(NCELL) ) GO TO 71 
9 CALL PATH! 10, IF) 

7 S = OUAOS( 10, IF ) 

GO TO 30 

71 CALL X ITP( ID) 

CALL P^THi 10, ITP) 

8 S = QUAOTP(ITP,IO) 
lOl = 10 

IF(FPATH.LT.2.0^S) GO TO 10 
81 T = SIGMA^(XTP-XD( ID) ) 

IF(T,LT,0,0) GO TO 82 

TALLY = 2.0/SQRT( USQG + TPHIOUO) I 

Y( ID) = Y( ID ) + TALLY 

N( ID) = N( ID) +2 

10=10+ IF IX( SIGMA ) 


44 


GO TO 81 

82 SIGMA = -SIGMA 

ID = lOL + IF IX ( SIGMA ) 

FPATH = FPATH - 2.0=s^S 
IF(XO.NE,XB( IBO) ) GO TO 33 
I8F = IBO 
RETURN 

33 I = IBG 
IBO = I8F 
IBF = I 
XF = XB( IBF) 

IF = NOEL<^( IRF-l ) ^ L 
IFt in.NE. IIP ) on TO 7 
S = OUAOTPI 10, IF ) 

GO TO 30 

10 IF(FPATH.GT.S) GO TO I? 

FPATH = S - FPATH 
CALL XICTPI ITP, 10) 

GO TO 41 

12 T = SIGMA^(XTP-XO( ID) ) 
IF(T.LT.O.O) GO TO 13 

TALLY = I.O/SQRT(USOG+TPHintIO) ) 
Y( ID) = Y( ID ) + TALLY 
N( ID) = N( ID ) + I 
ID = ID IF IX( SIGMA ) 

GO TO 12 

13 SIGMA = -SIGMA 

ID = ID 4^ IF IX(SIGMA) 

FPATH = FPATH - S 

I = IBO 

IBO = IBF 

IBF = I 

XF = XB{ IBF) 

IF = NDEL*{ IBF- 1) + 1 
CALL X ICTPI ITP, 10) 

GO TO 41 
END 
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tlBFTC STOSS 


SUBROUTINE STOSS 

COMMON /CM A IN/ NO , N 1 1 ALPH A , CONST ♦ NS t A I ( 20 ) ,VELI 1024 ), NK ♦ MFP ( 126 ) t 
« KFLAG, THETA, 0ISTB(64t64) tLFLAG 
COMMON /Cl TER/NI20) , Y( 20} , KNTR , VSQ f F P A TH , I BO, I BF , S I GMA , 1 0 , CRRNT 
COMMON/CMNPHI /NDELfOELXfXB ( 201 , TPHI 0(3 ,20) , I M I N ( 20 ) , PH I MI N ( 20 ) , 

♦ XMIN( 20) , TPHIOI 20) ,TPHIX( 1026) 

COMMON /CCELL/ NC ELL , XO , XF , I 0, I F , USQO 
COMMON /CXICTP/ IC,XC 
COMMON /CPATH/ EV 
DATA TWDP 1/6. 2B31853/ 

WSQ=U SQO+VSQ+TPHI X( 1C ) 

COSO ^ SIGMA^'SQRTr (USQU + TPHIXC IC) ) /WSQ) 

CALL RANO(R) 

IF(LFLAG.EQ.2) GO TO 3 
CnSl = 1.0-2.0*R 
GO TO 4 

3 I = IF IXC 64.0»R ) 4- 1 

IF( EV.GE.1.25) GO TO I 
J = IFIX( 16.0*EV)+1 
GO TO 6 

1 IF(EV.GE. 12.25) STOP 

J = IFIX(4.0«CEV-I.25) )+2l 
5 COSl ^ DI STB{ I , J ) 

4 CALL RANO(R) 

COSN = CO SO*COS 1-COS { TWOPI*R)*SQRTC ( 1. 0-C0S0*»2 ) ♦( 1.0-COS1*«2) ) 
USO = WSQ<^C0SN^«2 
USQO = USQ-TPHIX( IC ) 

VSQ = WSQ-USQ 
SIG = SIGMA 

SIGMA = SIGN( SIGMA, COSN) 

IF( SIG«SI GMA.GT.0.0) GO TO 2 
ID ^ 10 ^ IF IX( SIGMA) 

I = IBO 
IBO = IBF 
IBF = I 
XF = X8( IBF) 

IF = NDEL«{ I8F-II 4- I 

2 10 = IC 
XO = XC 
RFTURN 
END 











SIBFTC PATH 


SUBROUTINE PATH! II, 12) 

COMMON /CPATH/ EV 

COMMON /Cl TER/N( 20) , Y{ 20) ,KNTR, VSO , F PA TH , I 30 , I BF, SIGMA, IO,CRRNT 
COMMON /CMAIN/ N 0 , N 1 , ALPH A , CONS T , NS , A I I 20 ) , VE L ( 102 ^ , NK , MF P ( 1 26 ) , 
« KFLAG, THETA, 0ISTB(64, 64) tLFLAG 

COMMON /CMNPHI/N0EL,r)ELX,X8(20) , TPH I Q ( 3 , 20 ) , I M I N ( 2 0 ) , PH I M I N C 20 ) , 

* XMIN(20),TPHIQ(20),TPH1X( 1026) 

COMMON /CCELL/ NCFLL, XD,XF , 10, I F, USQO 
REAL MFP 
CALL RANDIR) 

J = IF IX( 1024.*R ) + l 
FPATH = ALPHA«VEL(J) 

FV = THETA*( USQO+VSO + 0, (TPHI X( II ) 4-TPHIXU2) ) ) 

IF( FV.LT.0,0) STOP 
IF(KFLAG.EQ. 1) RETURN 
IF(EV.GE.UO) GO TO I 
J = IFIX( 100.0<^EV)4-l 
GO TO 3 

1 IF( FV,GE*i3.25) GO TO 2 

J = IFIX( 2.0*(EV-U0) )+lOL 
GO TO 3 

2 J = 126 

3 FPATH = FPATH«MFP(J) 

RETURN 

END 










iCBFTC XIC 


SUBROUTINE X IC i 10, IF ) 

COMMON /CXICTP/ IC,XC 

COMMON /C I TER /N( 2 0) ,Y( 20) , KNTR , VSQ , F P ATH , I BO , I BF , $ I GMA , IO,CRRNT 
C0MMDN/CMNPHI/NDEL,DELX,X8 (20) , TPHIQ(3 ,20) , I M I N ( 20 ) , PHI MI N( 20 ) , 
* XMIN ( 20) , TPHID( 20) ,TPHI XC 1026) 

KNTR = KNTR ^ I 
F = FPATH 
II = in 
13 = IF 

1 12 = (I3-Il)/2 4^ II 

2 IFC I2.EQ.il) GO TO 4 
S = OUADS( II, 12) 

IF(S.LT.F) GO TO 3 
13 = 12 

GO TO 1 

3 F = F ^ S 
II = 12 
GO TO I 

4 XC -= UFLX^FLOAK Il-l) 

IC = I 1 

RETURN 

END 









$1BFTC XICTP 


SUBROUTINE X IC TP ( I TP , I 0 ) 

COM.^ON /CXJCTP/ IC,XC 

COMMON /CITER/N( 20) ,YI20) ,KNTR, VSQfPPATH, I RO, I BF, SIGMA, ID,CRRNT 
COMMON/CMNPHI/NDEL ,0ELX,X8 (70) , TPH I Q ( 3 ,20 ), I M I N ( 20 ), PH I M I N ( 20 ) , 
* XMINt 20) , TPHI0(20) ,TPHIX( 1026) 

KNTR = KNTR + 1 
F = FPATH 
K = 1 
II = I TP 
13 = 10 

1 12 = ( 13- I 1) /2 + II 

2 IF( 12.EQ.I 1) GO TO 7 
GO TO ( ^,4),K 

3 S = OUAOTPI 11,12) 

IF(S,EQ*0.0) GO TO 7 
GO TO 5 

^ S = OUAOS( 11,12) 

5 IF(S*LT.F) GO TO 6 
13 = 12 

GO TO 1 

6 F = F - S 
II = 12 

K = 2 
GO TO 1 

7 XC = DELX*FLOAT( I l-I) 

IC = 11 

RETURN 

END 















SIBFTC XITP 


SUBROUTINE XITP< 10) 

COMMON /CXITP/ ITP, XTP 

COMMON /CCELL/ NCELL t XO, XF , ZZZZZZ , I F , USQO 

COMMON/CMNPHl/NOELf OELXfXB ( 20 ) » TPH I Q ( 3 1 20 ) 1 1 MI N ( 20 ) ^ PH I MI N ( 20 ) , 
* XMINI 20) , TPHIDI 20) ,TPHIXU026) 

11 = 10 

12 = IMIN(NCELL) 

1 M=( I2-IU/2 
IF(M.EQ.O) GO TO 5 
I=Il+M 

TEST=USQO+TPHIX( I ) 

IF( TEST)2.4, 3 

2 £ 2=1 

GO TO 1 

3 £1=1 

GO TO I 

4 ITP=I 
GO TO 6 

5 £TP=Il 

6 XTP = OELX«FLOAT{ ITP-n 
RETURN 

END 



Figure 20. - Flow chart for subroutine QUAD, 
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SIBFTC QUAD 


I 


FUNCTION QUADIZZZZZZ) 

COMi'ION /CM A IN/ NO, Nl , ALPHA , CON ST, NS ,AI (ZO) ,VEL ( 1024 ) , NK , MF P( 126 » , 

* KFLAG,THETA,0ISTB(64,64) ,LFLAG 

COMMON /Cl TER/N (20) ,y ( 20) ,KNTR, VSO,FPATH, I BO, I BF,S IGMA, ID, CRRNT 
COMMON /CC6LL/ NCELL , XO , XF , 1 0 , I F ,USOO 

CQMM0N/CMNPHI/NDEL,DELX,XB(20) ,TPH I 0(3,20) , I MI N ( 20 ) , PH I MI N( 20 ) , 

* XM IN ( 20) , TPHID( 20) ,TPH1X( 1026) 

DIMENSION H( 2) 

DATA H/5.65555556E-1, 8. 88888880E-1/ 

F(Y) = SORT! 1.0FVSQ/(USQO + Y)) 

QUAD = H( 1 )<'F( TPHIQ( 1, NCELL) ) + H ( 2 ) *F ( TPH I Q ( 2 , NCELL ) ) + 

* H( 1)*F( TPHIQI3, NCELL) ) 

QUAD = ABS(QUAD/( 2.0*FLOAT(NS) ) ) 

RETURN 

END 
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Figure 21, - Flow chart for subroutine QUADTP. 












$IBFTC QUAOTP 


FUNCTION QUADTPI IIP* IF) 

C OPFN ENDED INTEGRATION 

COMMON /CITER/N(20),Y( 20) tKNTRf VSQ*FPATH, I BOt 1 BF ,S I GMA , ID t CRRNT 
COMMON /CCELL/ NCELL , XQ , XF , I Ot 2 ZZZZ Z , USQO 

COMMON/CMNPHI/NOELfQELXtXB (20),TPH1Q(B ,20) , I MI N( 20 ) , PH I MINI 20 ) , 
♦ XMINI 20),TPHIO(20),TPHIX( 1026) 

DIMENSION HI 3) 

DATA H/4. 84974226,- 3. 9 191 8350^ 2 *4/ 

GtY) = SQRTI 1.0<-VSQ/(USQ0+Y) ) 

QUAOTP = 0 
I = lABSI ITP-IF) 

IFd.LT.3) RETURN 
J = 1/20 

IF(J.EQ*0) J = 1 

IF( ITP^GT.IF ) J = - J 

INI = ITP + J 

IN2 = INI + J 

IN3 = IN2 + J 

QUAOTP = HI 1 }*G( TPHIXI INI ) ) + H ( 2 ) *G ( TPH I X ( IN2 ) ) + 

^ HI 5)*GITPHIXI IN3) ) 

QUAOTP = A8SI FtOATIJ )*OELX«QUADTP) 

QUAOTP = QUAOTP + QUADStIN3,IF) 

RETURN 

END 



( quads) 



Figure 22. - Flow chart for subroutine QUADS. 
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$IBFTC QUADS 


FUNCTIUN QUADS(If),IF) 

COMMON /CITER/N< 2 0),Y( 20) , KNTR , VSQ » F P ATH ♦ I BO , I 8F # S I GMA , IO,CRRNT 
COMMON /CCELL/ NCELL # XO, XF f ZZZZ2Z C2 ) ,USOO 

COMMON/CMNPHI/NDELtOELXtXB t 20) , TPHIQt 3 ,20 ) , I MI N ( 20 ) , PH I M I N C 20 ) t 
^ XMIN{ 20) ♦ TPHIO< 20) fTPHlXC 1026) 

C SIMPSON’S RULE 

F(Y) = SQKT( l*0+VSQ/(USQn^Y) ) 

QUADS = 0 

IF< lO.EQ* IF) RETURN 
FI = MINO( 10, IF) 

12 = MAXOF in, IF ) 

13 = 12-2 

8 J < I2-I I )/10 
IF(J.EQ.O) J=1 

OELXJ = FLOAT! J )^DFLX/3,0 

9 DO 10 I=Il, 13,2 
INI = ( I- I l)*J+I 1 
IN2 = INl+J 

IN3 = IN2+J 

IF( I .GT.I 3) GO TO 12 

IF( IN3.GT. 12 ) GO TO ll 

10 QUADS = QUADS + DEL XJ» ( F ( TPH I X (I N1 )) +4. 0*F ( TPH I X ( I N2 )) + 

♦F( TPHIXC IN3) ) ) 

QUADS = ABS(QUAOS) 

RETURN 

11 II = INI 
GO TO 8 

12 QUADS = ABSFQUADS ^-OELX* (F ( TPHI X( 12-1 ) ) +F ( TPHI X( 12) ) )/2.0) 

RETURN 

END 







$IBFTC CHE8Y 


SUBROUTINE CHEBY(MODF) 

COMMON /CCHEBY/ B(?0)tXD120) 

COMMON /CM A IN/ NOf N , ALPHA ,C ONS T ,NS , A I ( 20) tVEL(1024> ,NK , MF P( 1 26 ) » 

* KFLAG,ZZZZZ,0ISTB(64t64) ,LFLAG 

COMMON /Cl TER/ZZZI 20) , Y(20) ,KNTR,VSO»FPATH ,I BOt I 8F , S IGMA» lOtCRRNT 
DIMENSION T(20,20),C(20) 

DOUBLE PRECISIUN PI, THETA 
DATA P 1/3. 141 59265 35 897932/ 

GO TO <10, 20), MODE 
10 EN=N-1 
NMl=EN 
RN=1./EN 
TRN ^ 2.«KN 
I = 0 

X0( I ) = - l.0E^30 
DO 12 I=l,N 
A=I-l 

THETA= A*RN*P I 

12 XO(I) = .5<=( I.-DCQS(THETA) ) 

XD(N+1 ) = 1. OE+30 
DO 15 1=1, N 

DO 14 J=1,N 
R^J-1 
S=N-I 

THETA = R^S^^RN^PI 

14 T( I,J )= DCOSl THETA) 

15 CONTINUE 
RETURN 

20 CONTINUE 

DO 26 K=1,N 

C(K ) = .5^( Y< 1)«T( 1,K) + Y(NH^T(N,K)) 

DO 24 1=2, NMl 

24 C(K)= CiK)^ Y( T(I,K) 

26 8(K) = TRN*C(K ) 

B(N) ^ B(N)/2.0 

RETURN 

END 








SIBFTC PHI 


FUNCTION PHKYI 

COMMON /CHAIN/ NO»N, ALPHA , CONST ^NS , A I (20) ,VEL( 1024) tNKf MFP( 126 ) t 
♦ KFLAGtTHETAf DISTB(64t64) fLFLAG 

dimension 8B(22) 

X=2.*Y-1, 

B8(N+2) = 0 
30(N+l) = 0 
C=2** X 
DO 10 J=lfN 
K = N + l-J 

10 B8(K) = C^B8(K + l)-B8(K+2) AI(K) 

PHI = (8B( 1)-8B(3) )/2.0 

RETURN 

END 








$18FTC OPHI 


FUNCTION DPHKYI 
COMMON /CCLN2S/ DAI20> 

COMMON /CMAIN/ NO, N t ALPHA .CONST ,NS t A I( 20) ,VEt (1024) ,NK,MFP( 126), 
♦ KFLAG, THETA, 0ISTB(64, 64) ,LFLAG 
DIMENSION BB(2?) 

X=2.*Y-l. 

B8(N+2) = 0 
B8(N+l) = 0 
C=2.* X 
DO 10 J=1,N 
K = N4-1-J 

10 BB(K) = C*BB(K+l )-8B(K*2) + OA(K) 

DPHI = (BB(1)-B8( 3) )/2.0 

RETURN 

END 
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$IBFTC 


DFNS 


FUNCTION OENS(Y) 

COMMON /CCHE8Y/ B(20)tX0(20) 

COMMON /CMAIN/ NO, N, ALPHA , CONST, NS, A I ( 20) , VFL( 1024) , NK , MF P ( 126 ) , 
♦ KFLAG, THETA, 0ISTB(64, 64) ,LFLAG 
DIMENSION BB(22) 

X = 2*x^Y-l. 

B8(N+2) = 0 
8B(N+1) = 0 
C=2,« X 
on 10 j=i,N 

K = N4-1-J 

10 BB(K) == C#BBC K+l)-BB(K+2) + B(K) 

DENS ^ (BBU)-B8(3))/2.0 

RETURN 

END 
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SIBFTC OISCRl 


SUBROUTINE D I SCR I ( M ,N f A , ME AN A, STOA ) 

REAL MEANA 

DOUBLE PRECISION FM,MEAN1 , SUMl » SUM2 
DIMENSION A( 20,20) ,MEANA ( 20) , STOA(20) 

FM = M 
FMl = M-l 
00 2 J=l,N 
SUMl = 0.000 
SUM2 = O.ODO 
DO I 1=1, M 

SUMl = SUMl + A( I, J) 

1 SUM2 = SUM2 + AU , J)<=A( I ,J) 

MEANl = SUMl/FM 

STOA(J) = SQRT(ABS( (SUM2/FM-MEAN1=^^MEANI)/FM1) } 

2 MEANA(J) = MEANl 
RETURN 

END 
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Figure 28. - Flow chart for subroutine DISCR2. 
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SIBFTC DISCR2 


SUBROUTINE 0 ISCK 2 ( K , A , AMEAN , A06 V) 

DIMENSION AC20) 

DOUBLE PRECISION FKR , iME AN 1 f SUMl , SUM2 
FKR = l*ODO/DBLE(FLOAT(K) ) 

SUMl = 0. 

SUM2 = 0. 

DO 13 J=1,K 

SUMl -= SUMl + A(J ) 

13 SUM2 = SUM2 + AC J ) *A ( J) 

MEANl =SUM1’«‘ FKR 

ADEV = SQRTC ABSCC SUM2*FKR-MEAN1*MEANI) /FLOATCK-l) ) )/ 

AMEAN = MEANl 

RETURN 

END 










o r> n o n n 


ilBFTC PLOTF 


SUBROUTINE PLO TF ( Nf XO , XF ,F OFX ) 

SUBROUTINE FOR SINGLE-PAGE PLOTTING OF FUNCTION FOFXIXI FROM XO TO XF 
USING N POINTS (N OOO.ANO.LE. lOl) 

CALLS SUBROUTINE PLOTYX 


DIMENSION X( 101) ,Y ( 101 ) 

YMAX = -I.FIO 
YMIN = I.EIO 
DEL=(XF-XO)/FLOAT(N-l) 

DO 10 1=1, N 

X(I) = OEL*FLOATI I-l) ♦ XO 

Y( 1 ) = FOFXIXI I J I 

IFl YI I I.GT.YMAX) YMAX = YII) 

10 IFIYl I I.LT.YMIN) YMIN = YII) 

CALL PLOTYXI N , Y, X, YMAX , YMIN , XO, XF » 

RETURN 

END 
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o o o o n o o 


$IBFTC PLOTYX 


SUBROUTINE PLU TYX( N , Y, X , YM AX , YMI N ,XO , XF ) 

SUBROUTINE FOR N-POINT, SINGLE-PAGE PLOTTING OF ARRAYS X(NJ,Y(N) 

N MUST BE OOO.ANO.LE. lOl 

CALLS SUBROUTINE SORTYX ANO SCALEY 

DIMENSION Y( N) f X(N ),P ( m 
DATA P/ll*0./ 

CALL SORTYXIN , Y,X) 

CALL SCALEYIN, Y,YMAX,YMIN,P(6) ,P17> ,P ( 8) ,NSCALE ) 

Pin = N 

NN = ALOGlO( P( l)-l, ) 

P(9) =0 

P( 10) = l.E6*X0 

PI 11) = IXF-XO) * I.E4 

CALL PLOTXYl Y,X,ll8tP) 

WRITE(6,10) NSCALE 

10 F0RMAT(2HPLf20X, 35HSCALE FAC TOR F OR ■ ORDI NATES IS 10** ( , I 3, I H) ) 
RETURN 
END 











SIBFTC SORTYX 


SUBROUTINE SORTYXC Nt Y # X I 
DIMENSION Y(N)fX(N) 

M=N-l 

DO 10 I = l,M 
K = I + l 

00 10 J ^ K,N 
A = Y( 1) 

B = x( n 

IFIA.GF.YI J) I GO TO 10 
Y( II = Y( J I 
XI n = XIJ) 

YU) = A 
XIJ) = B 
10 CONTINUE 
RETURN 
END 



( scaley) 



Figure 32. - Flow chart for subroutine SCALEY. 





I 


SIBFTC SCALEY 

SUBROUTINE SCALEY ( N, Y , YMAX , YMI N ,KSY ,F Y .DY, NSCALE ) 
DIMENSION YIN) tNOY (20) ,NSYI20) 

DATA N0Y/2f 4,625,2*1, 125,1429,3*2.3*25,4*3333 t 3*5/ 
DATA NSY/ 2*4, 2, 2*5, 3, 2, 3*5, 3*4, 4*2, 3*5/ 

REAL KSY 
2MAX=ABS( YMAX) 

ZMIN=ABS(YMIN) 

NSCALE = ALOGIOIAMAXUZMAX.ZMIN)) 

A =10.**(-NSCALE) 

NPLUS = YMAX * A 

NMINUS = YMIN * A 

IF! YMAX .GT.O.) NPLUS=NPLUS+1 

IFINMINUS.LT.O ) NMI NUS=NMI NUS- 1 

J = NPLUS - NMINUS 

KSY = NSYIJ) 

DY =-NOY(J) 

FY = NPLUS *10*»(6-NSY(J) ) 

00 10 I = 1,N 
10 Y( I ) = A*YI 1 ) 

RETURN 

END 
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Auxiliary FORTRAN IV Program Descriptions 

The programs given in this section are necessary to generate the binary formatted 
data decks discussed in the sections Tabulated Distributions and Preparation of Input 
Tables. For an example of input to and printed output from these programs, see appen- 
dix F, The ranges of the independent variables used in programs MFP, ARGON, and 
ARGINV were chosen specifically for the gas used in the sample problem. If the user 
wishes to change these ranges, for a different gas, he must also change the subroutines 
PATH and STOSS of ENEC accordingly. In other words, ENEC expects the tables to be 
given for specific ranges. 


CVEL 


Purpose: 

To construct and punch on cards a table of the functional values of -In R for 1024 
equally spaced values of R over the range from 0 to 1 
Method: 

The table is punched on cards by subroutine BCDUMP. 

Program called: 

BCDUMP (appendix E) 

FORTRAN listing: 

SIBFTC CVEL 

DIMENSION VEH 1024) 

DELX = 1.0/1024.0 

DO I 1=1, 1024 

X = DELX*( FLOAT! I-H+0.5) 

1 VELin = -ALOG(X) 

CALL BCDUMP! VFL! 1 ) ,VFL! 10241 ) 

DO 2 L=l,4 
WRITE !6,201) 

DO 2 1=1,256,8 

11 = H-1L-1)*256 

12 = 11+7 

WRITE 16,202) 1 1 , ( VEL ! J ) , J = I 1, I 2) 

2 CONTINUE 

201 FORMAT ! IHl, /// , lOX , IHI , 44X , 3H VEL , //) 

202 FORMAT !IH , 5X, I 5 , 8F 11 . 4) 

STOP 

END 
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MFP 


Purpose: 

To construct and punch on cards a table of A(E) = 1.01 /ct(E) for 100 values of E 
equally spaced over the range from 0 to 1 and 25 values equally spaced over the 
range from 1 to 13. 5 
Method: 

1/2 

The values of X(E) are obtained by curve fitting known values of a(E ' ) with a cubic 
spline (appendix B) and interpolating. The last value in the table is set to 0. 012 and 
is used by ENEC for energies greater than 13. 5 electron volts. The table is then 
punched on cards by subroutine BCDUMP. 

Programs called: 

BCDUMP (appendix E) 

SPLINE and related subprograms (appendix B) 

FORTRAN listing: 

»IBFTC MFP 

C COMPUTES MFP TABLE FROM SPLINE FIT OF SIGMA S 
REAL MFP 

DIMENSION SQTEV(21),SIGMASI21) ,MFP( 126 J ,E f 125 ) 

READ (5,10U (SQTEVID.SIGMASII ),I = l»21) 

WRITE ( 6 , 201 J ( SQTEVI I ) .SIGMASm ,1=1,21) 

CALL SPLINE! SQTE V, SIGMAS , 2 1 , 0. 0,0. 0 ,2) 

00 2 1=1,100 

2 Ed) = .005+.01’tFL0ATI I-l) 

DO 3 1=101,125 

3 Ed) = 1.25^-.5*FLOATII-101) 

DO 4 1=1,125 

4MFPII) = l.Ol/FI SORTIE I I ) ) ) 

MFPI126) = .012 

CALL BCOUMPIMFPI 1) ,MFP(126) ) 

WRITE (6,202) d , E d ) , MF PI I ) , I = 1 , 1 2 5 ) 

I = 126 

WRITE 16,203) I,MFP(126) 

STOP 

101 FORMAT I 2E10.0) 

201 FORMAT I IH 1 , /// , 1 IX, 8HSQRT ( E V) ,10X, 7HS IGMA-S, // , I F20.7 , F15 . 2 ) ) 

202 FORMAT (IH 1, /// , 20X, 2HE V , I 3X, 3HMFP , // , 1 1 4, F20. 3 ,F20 .8 ) ) 

203 FORMAT d 4, 20X , F20. 8 ) 

END 
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Figure 33. - Flow chart for subroutine MFP. 


ARGON, G, SIMPS 

Purpose: 

To compute values on the surface defined by equation (11) for input to program 
ARGINV 
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Method: 

These surface values are computed by surface fitting (appendix B) known values of the 
differential cross sections a(0,E) for the gas under consideration and performing 
Simpson’s Rule integration. 

Remarks: 

ARGON uses two function subprograms G and SIMPS to perform the Simpson’s Rule 
integration 
Programs called: 

SIMPS 


SPLIN2 and related subprograms (appendix B) 

FORTRAN listings: 

$IBFTC ARGON 

DIMENSION X( 20) »Y( 20) ,Z (20,20) fXl(20) , Zl (20,20) ,SIGS< 20), S( 20,20) 
COMMON /FUN/ YY 
EXTERNAL G 

READ ( 5, lOl) ( Xl( I) ,1 = 1, 13) 

READ (5,102) <Y(J),(Zl(I,J),I=l,l3) ,J = l,20) 

00 I 1=1, 13 
K = 14-1 

X(I) = C0S(Xl(K)/57*2957795) 

00 I J=l,20 

1 Z( I,J ) = ZKK, J) 

CALL SPLIN2CX,Y,Z, 13,20) 

00 2 J=l,20 
YY = Y( J ) 

2 SIGS(J) = SIMPS(-l*0,1.0,G) 

00 3 1=1,11 

X(I) = -l.0i-0.2<'FLnAT( I-l) 

DO 3 J =1,20 
YY = Y(J ) 

3 S(I,J) = SIMPS(-1.0,X( I ) ,G)/SIGS( J) 

WRITE (6,203) ( X ( I ) , 1 =1 , 1 1 ) 

on 5 J=l,20 

5 WRTTF (6,204) Y( J ) , ( S ( I , J ) , 1 =l , 1 1 ) , SI G S ( J ) 

STOP 

101 FORMAT ( 5X ,I3E5.0) 

102 FORMAT ( 14E5.0) 

203 FORMAT ( IH 1 , // /, 60 X, 5H S I GM A , /, 5X , I 3HEV /COS ( THE T A ) , 1 1F8 • 1 , 2X , 7HS IGM 
*A-S,//) 

204 FORMAT I 3X , F 5 • 2, 1 OX , I IF 8. 4,F 8* 3 ) 

END 

$I8FTC G 

FUNCTION G(X) 

COMMON /FUN/ YY 
G = F(X,YY) 

RETURN 

END 
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SIBFTC SIMPS 


FUNCTION SIMPS(XMIN,XMAX,FI 

OELX = IXMAX-XMIN )/128*0 

SUM = 0.0 

DO I 1=1,128,2 

XI = XMIN+DELX*FLnAT( I-l) 

X2 = Xl+OELX 
X3 = X2+DELX 

I SUM = SUM4^F(XU+4.0<'F(X2)+F(X3) 
SIMPS = DELX*SUM/3.0 
RETURN 
END 
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Flow charts: 


C ARGON 0 


Read input 
data 







Figure 34. > Flow chart for subroutine ARGON. 


Figure 35. - Flow chart for subroutine SIMPS. 


ARGINV 


Purpose : 

To solve equation (11) with R being taken at 64 equally spaced points ov .r the range 
from 0 to 1 and 20 and 44 values of E equally spaced over the ranges 0 to 1.25 and 
1.25 to 12.25, respectively, and to punch this data on cards 
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I II III 


Method: 

The equation is solved by Newton-Raphson iteration of functional values obtained from 
a two-dimensional surface fit (appendix B) of the output from ARGON. Convergence 
of the Newton-Raphson iteration is assumed when the difference between two consecu- 
tive iterations is less than 0. 0001. The constructed table is punched on cards by 
BCDUMP. 

Programs called: 

BCDUMP (appendix E) 

SPLIN2 and related subprograms (appendix B) 

FORTRAN listing: 

$IBFTC ARGINV 

C FINOS THE INVERSE FOR THE ARGON SURFACE 

DIMENSION X(20)tYl20)»Z<20»20> ,R(64),E{64) ,ARG(64,64) 

READ (5,101) tX( I ),I = l,ll) 

READ (5,102) (Y(J),(Z(I,J),T=1,11) ,J=1,20) 

WRITE (6,201) 

CALL SPL IN2I X,Y,Z,ll,20) 

DELR = 1.0/64.0 
R( 1) = DELR/2.0 
DO 1 1=2,64 

1 R( I ) = R( I-D+OELR 
DELE = 1.25/20.0 
Ed) = OELE/2.0 

00 2 J=2,20 

2 E( J ) = E( J-D+OELE 

DELE = ( 12.25-1.25)/44,0 
E(21) = 1.25+DELE/2.0 
DO 7 J=22,64 

7 E( J ) = E( J-IH-OELE 
A = -0.8 

ERROR = l.OE-4 
00 4 1=1,64 
00 4 J = l,64 
EE = E(J ) 

00 3 K=l,60 

OELA = (F(A,EE)-R( I))/FX(A,FE) 

IF(ABS(OELA) .LT. ERROR) GO TO 4 

5 IF(ABS(A-0ELA).LT.1.0) GO TO 6 
DELA = OELA/2.0 

GO TO 5 

6 A = A-OELA 

3 CONTINUE 
WRITE (6,209) 

WRITE (6,210) I,.',Rn ) ,F( J),A,OELA 

4 ARG( 1, J) = A 

CALL BCOUMP(ARG( 1,1) ,ARG(64,64) ) 

WRITE (6,202) 

DO 8 1=1,64,8 

11 = I 

12 = 11+7 

8 WRITE (6,203) 1 1 , ( R( K ) , K=I I , I 2 ) 

WRITE (6,204) 

00 9 J=l,64,8 
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■ I III 


IIMII ■■ II 


HI HI II HIM ■III II II I I II I 






i 

I 

I 

t 


J1 = J 
J2 = Jl+7 

9 WRITE J6,205) J1 , ( E ( K ) , K=J I , J2 ) 

DO 10 L=l, 16 
WRITE (6,2061 

00 10 M=l,4 

1 = 4=^(L-1)+M 
WRITE (6,207) 

DO 10 J=l,64,8 
J1 = J 

J2 = Jl + 7 

WRITE (6,208 ) I, J1,(ARG(I ,K) ,K=J1, J2) 

10 CONTINUE 
STOP 

101 FORMAT (6X,11E6.0) 

102 FORMAT ( 12E6.0) 

201 FORMAT ( IHl) 

202 FORMAT ( IH 1 , /// , 5X , IHI , 50X , IHR , //) 

203 FORMAT ( IH , I 5 , 5X , 8F II . 7 ) 

204 FORMAT ( IHL , lOX , IH J ,45 X , 2HE V , / / ) 

205 FORMAT ( IH , 5X , I 5 , 8F 1 1 . 7 ) 

206 FORMAT ( IH 1 , / / , 5X , IH I , 4X , IH J ) 

207 FORMAT ( IHK) 

208 FORMAT ( IH ,2I5,8F10,4) 

209 FORMAT (54HKTHE FOLLOWLNG DID NOT CONVERGE AFTER SIXTY ITERATIONS) 

210 FORMAT (3HKI = ,I4,5X, 2H J = , I 4 , 5 X , 2HR = , E 1 7. 8 , 5X , 2H E= , E 17 . 8 , 5X , 2HA= , E 1 
♦ 7.8,5X,5HOELA = ,E 17.8) 

END 
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Flow chart 
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APPENDIX A 


SYMBOLS 


A 

c 

£ 

e 

f(u,x) 

f(u,V,x) 

AAA 

f(u,v,w) 

Gx(t) 

g(«,V) 

g(u,x) 

g(u,V,x) 

J 

Jo 

k 

L 

I 

M 


normalization factor 
space charge parameter 
electron kinetic energy 
electron charge 

marginal probability function of f(u,V,x) 
probability function 

Maxwellian velocity distribution function 

cumulative distribution function, P(X < t) 

nondimensionalized flux probability function 

marginal probability function of g(u, V,x) 

probability function 

current to collector 

emission current 

Boltzmann constant 

electrode separation 

path length along electron trajectory 

path length along electron trajectory to collision 

maximum value of dimensionless electron density distribution 


Nc 

No 

n(x) 

^g 

Po 

P(0 < 0) 
R 


electron mass 

electron histories that reach collector 
total electron histories for one iteration 

A , 

dimensionless electron density distribution, p(x)/p(0) 
scattering gas pressure 

273 

reduced pressure, P — ~ 

o ^ 

s 

probability that random variable 9 is less than or equal to some 0 
random number uniformly distributed over range from 0 to 1 





T 

T 

u 


g 


u,v,w 

V 


random number uniformly distributed over range from 0 to 1 associated with a 
random variable X 

emitter temperature 

temperature of scattering gas 

dimensionless velocity component in x-direction, u(2kT/m)^'^^ 
velocity components 


dimensionless velocity component transverse to u, |(2kT/m)(v^ - w^)J 


1/2 


X dimensionless spatial variable normal to electrodes, x/L 

A 

X spatial variable normal to electrodes 

Xp point of collision 

Xj^ Chebyshev abscissas for curve fit 

Xq location of last electron ’’event’* (cell boundary or collision) 

€ permittivity of free space 

6 polar angle and polar angle in laboratory system after collision, (eq. (12)) 

6^ polar angle in laboratory system at point of (before) collision, (eq. (12)) 

6' polar angle in center-of-mass system after scattering 

X(E) energy -dependent mean free path 

X(x) mean free path 

V potential distribution, (eq. (1)) 

v(x) computed potential distribution 

p(x) electron density distribution 

p(O)'*’ electron density of emitted electrons 

o (E) total scattering cross section 

a(0,E) differential scattering cross section 

(p azimuthal angle of incidence 

(p' azimuthal angle in center-of-mass system after scattering, (eq. (12)) 

<p{x) dimensionless potential distribution, ev(x)/kT 

Superscript: 

(—) average 
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APPENDIX B 


SPLINE CURVE AND SURFACE FITS 

It is the purpose of this appendix to describe the subroutines needed to perform one- 
and two-dimensional interpolatory spline curve fits. The subprograms presented herein 
are used by the programs that prepare the input tables for ENEC . 


ONE-DIMENSIONAL SPLINE CURVE FIT 

The specific functions used in the spline curve fit are piecewise polynomials of the 
third degree with hiatching first and second derivatives at the data points. These func- 
tions yield excellent approximations to the curve being fit as well as to the first deriva- 
tives of the curve. 


Method 

For a given set of n distinct data points x. in increasing order (x. < the 

corresponding functional values f(xj^), and either f*(xj) and f'(Xj^) or f"(Xj) and f”(Xj^), 
a cubic polynomial may be fit between each two data points subject to the following con- 
straints : 


g;’(Xi^i) = g{;i(Xi^i) 


(Bl) 


where g. denotes the cubic between the data points x^ and x.^j, and the primes denote 
differentiation with respect to x. From the Hermite interpolation formula (ref. 9) 


gj(x) = h.(x)f(x^) + hj^i(x)f(x.^j) + h.(x)f(x.) + hj^i(x)f’(x.^j) I 


Xj < X < x^^j i = 1, 2, . . . , n J 


(B2) 


and from the constraints given in equation (Bl), the following formula is obtained: 
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3 


(Xi^l - x.)f’(x..i) + 2(x.^i - x._j)f(Xj) + (X. - x._i)f’(Xi^l) = 


(x .^1 - x.)(x. - x._i) 


X |(X. - x._i)2f(x.^l) + [(x .^1 - x.)2 . (X. - x._i)2] 

X f(x.) - (x.^j - x/f(Xj_i)j i = 2, 3, . . . , n - 1 (B3) 

Given values Vj and Vg for f’(Xj) and f'CXj^), the following equations may be written: 

t’(Xj) = Vj 

f'(x„)=V2 ^ 

If the values for f”(xj) and f”(Xj^) are given, equation (B2) may be used to obtain 

3 


(B4) 


2f(Xj) ^ fCxj) = ^ ' [Kxj) - f(x,)] - i(X2 - XjjVj 


f(x„.j) + 2C(x„) = 


(X2 - X,) 

^ [«x„) - t(x„.i)] + i(X|j - x„.j)V^ 


> (B5) 


(=‘n - *n-l> 


where Vg and are the given values of f”(Xj) and f”(Xj^), respectively. 

The set of equations defined by equations (B3) and (B4), or the set defined by equa- 
tions (B3) and (B5), form a system of n linear equations that may be solved for f’(x.). 
This system may be written in matrix notation where the matrix to be inverted in solving 
the system is tridiagonal 
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where 


F’ = [f’(xj) f’(x2) . . . f’(xj] 


and 


^ " [Pl ' ‘ • ^n] 


Ai = xj^i - Xj 


Bi = 2(Xi^i - x._i) 


Ci = Xi - X. _i 


- ^+1 - 


i=2, 3, . . . ,n-l 


By specifying f’(xj) = Vj and f’(Xjj) = Vg, 


An = 0 


Bi = Bn = 1 


Cj = 0 


Di = Vi 


= Vo 
n ^ 


By specifying f”(xj) = Vg and = V^, 




Bi = Bn = 2 


Ci = l 
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3 


D 


n 




(X2 - Xj 

_3 

- *n-l> 


- [f(x2> - f(xi>] - i(x2 - Xj)V3 

[f(x„) -f(Xn_i)]+l(Xn-Xn_i)V4 


This tridiagonal system may be solved by the following algorithm according to Peaceman 
and Rachford (ref. 10) . 

Let 


Wi = Bi 


Wi = Bi - A.b._i 


i=2, 3, 


n 



i = 1, 2, . . . , n - 1 


d 


1 



di = 


Pj - ^^i-1 

W, 


i — 2y 3j • • • ^n 


then 


f’(x„) = d„ 

f’(xj) = d. - bjf’Cxj^j) i = 1, 2, . . . , n - 1 

Now that all the f’(x^) have been determined, the desired set of cubic interpolation equa- 
tions ^(x) may be obtained by returning to the Hermite interpolation formula equa- 
tion (B2) and deriving equation (B7) : 
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g.(x) = (x - X.)* 


2f(xp - 2f(x.^j) + (x.^j - x.)f'(xp + (x.^j - Xj)f'(x.^j) 


^+1 - 


+ (x - xj' 


-3f(x.) H- 3f(x.^^) - 2(x.^i - x.)f»(xp - (x .^1 - x.)f»(x.^^)~ 

(Xi^l - xf 


+ (x - x.)f’(xp + f(xp X. < X < x .^1 (B7) 

A set of second-order polynomial interpolation equations for f'(x) may be obtained by 
simply differentiating equation (B7) with respect to x. 


FORTRAN IV Subprograms for SPLINE Curve Fit 


Subroutine SPLINE (XX, Y,NN, Bl, BN, J): 

XX Array of independent variable in increasing order 

Y Array of functional values of curve to be fit; order must correspond to XX 
array 

NN Number of values in XX array; NN < 100 
Bl Boundary condition at XX(1) 

BN Boundary condition at XX(NN) 

J Boundary conditions specified are first derivatives; J = 1 

J Boundary conditions specified are second derivatives; J = 2 

Purpose: 

To construct and solve the tridiagonal system of equations (eqs, (B6)) and to compute 
the coefficients of the interpolation equations gj^(x). These coefficients are stored in 
COMMON/SPLN/. 

Labeled COMMON: 

/SPLN/A(100), B(IOO), C(IOO), D(100),X(100), N 

A,B,C,D Arrays of coefficients of zero, first, second, and third degree terms of 
interpolation equations (eq. (B7)) 


N 


Array of independent variable corresponding to user’s XX array 
Number of values in X array corresponding to user’s NN 
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FORTRAN listing: 

S18FTC SPLINE 

SUBROUTINE SPL INE ( XX , Y ,NN , B 1 , B N t J) 

DIMENSION XX{ 1),Y( 1) tOYI 1001 

COMMON /SPLN/ A ( 1 00 ) , B < 1 00 ) tC ( 100) , 0 C 1 00 I , XIIOO J , N 
N = NN 
DO 7 1=1, N 

7 X( I ) = XXI I) 

Nl = N-1 

GO TO (4,5),J 

4 Bt 1 ) = 1,0 
C( 1 ) = 0.0 
D( 1 ) =81 
A(N) = 0.0 
B(N ) = 1.0 
D(N) = BN 
GO TO 6 

5 BC 1) = 2.0 

C{ 1) = 1.0 

Oil) = 3.0^(Y(2)-Y( I))/(X<2)-X{1))-0.5*(XC2)~X(1))#'R1 

A(N) = 1.0 

B(N) = 2.0 

0(N) = 3.0=«^l Y{ N)-Y(Nl) ) /( Xf N)-X<N1 ) X(N)-X(Nin *8N 

6 DO 1 1=2, Nl 

Ad) = Xd+1 )-X( I ) 

B( I ) = 2.0«(Xt I+l)-X( dl) ) 

C( I) = X( I )“X< I-l) 

1 0(1) = 3.0*(Y( I+I)’^=C(I)’«^*2^-Y(n*CAn)’«'^2-Cm**2)-Yn-l)*A(I)**2)/ 
♦ (C( I)»A{ I ) ) 

D( 1) =0(1 )/8( 1) 

DO 2 1=2, N 

B( I) = B( I )-A( (--1)/8(I-1) 

2 0(1) = {0( i)-A( i)*on-i))/B(n 
DY ( N ) = 0 ( N ) 

DU 3 I = 1 , N 1 
K = N-I 

3 OY(K) = D( K)-C (K)«OY(K+U /B (K) 

A( I ) = X( 2 )-XI 1) 

00 8 1=1, Nl 

D(I) = ( ?.0*Y( I )-2.0*Y( I + I ) + A< I )^nY(I)>A( I )^0Y( I+I) )/A( n^*3 
C( I) = (-3.0<^Y( I ) + 3.0=i'Y( 1+ l)-2. 0*A( I )*DY( d-A( I )«0Y( I + l) )/AII )«^2 
B( I ) = DY( I ) 

8 A( 1 ) = Y( I ) 

RETURN 

END 

FUNCTION F(X1): 

XI Independent variable XX(1) < XI < XX(NN) 

Purpose: 

To apply the coefficients in COMMON/SPLN/ and to determine the interpolated value 
of the curve fit at XI by using equation (B7). 

Labeled COMMON: 

/SPLN/ 
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FORTRAN listing: 

SIBFTC F 

FUNCTION FCXl) 

COMMON /SPLN/ A ( I 00 ) , B ( 1 00 ) ,C ( 100) ,0 ( I 00) , X ( 1 00 I ,N 
1F(XI ,LT,Xl in STOP 
on 1 I=2,N 
J = I-l 

1 IF( XI ,LE,X( I ) > GO TO 2 
STOP 

22= Xl-X( J ) 

F = A{J)+Z*{B(J)+Z*(C(J)+Z*0(J))) 

RETURN 

END 

FUNCTION DF(X1) : 

XI Independent variable XX(1) < XI < XX(NN) 

Purpose: 

To apply the coefficients in COMMON/SPLN/ and to determine the interpolated value 
of the first derivative of the curve fit at XI by using the first derivative of equa- 
tion (B7). 

Labeled COMMON: 

/SPLN/ 

FORTRAN listing: 

$I8FTC OF 

FUNCTION DF(Xl) 

COMMON /SPLN/ At 1 00 ) , 3 { 1 00 ) ,C ( 1 00 ) , 0 ( I 00 > , X (1 00 ) t N 
IF( XI .LT.Xt 1) ) STOP 
DO 1 1=2, N 
J = I-l 

1 IFtXl .LE.Xt I) ) GO TO 2 
STOP 

2 Z = X 1-Xt J ) 

OF = B(J)+Z*(2.0*C(J)+3.0*?*0( J») 

RETURN 

END 

Program use: 

The user must call subroutine SPLINE, with the proper arguments, only once for the 
curve to be fit. After this call, the user has available to him labeled COMMON 
/SPLN/ and may use both functions F and DF. The user must be sure that F and 
DF are only used with an argument in the range of the independent variable because 
of stops built into the functions. This curve fit should not be used for extrapolation. 


TWO-DIMENSIONAL SPLINE SURFACE FIT 


The two-dimensional spline is completely analogous to the one-dimensional case in 
that the surface to be fit, defined on a rectangular grid with boundary conditions along the 
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edges, is redefined in terms of one-dimensional spline curve fits along the grid lines 
parallel to the coordinate axes. Two-dimensional interpolation is achieved by multiple 
one-dimensional interpolations of the function and its derivative. 


Method 


The method of solving the tridiagonal system of equations and obtaining the interpola- 
tion equations is identical to the one -dimensional case and will not be repeated; however, 
the boundary conditions that must be supplied need clarification. The surface to be fit 
must be defined by a set of n x m functional values on a n x m rectangular grid. Also, 
boundary conditions consisting of first or second partials (of the function being fit) must 
be supplied for the four boundaries. For example, if f(x,y) is the surface to be fit, it 
must be defined by functional values f(x.,y.) on the rectangular grid x^ x y^ for 
i = 1, 2, . . . , n and j = 1, 2, . . . , m, and by one of the following sets of boundary 
conditions : 





^9y/x=x. 

y=y. 


m 


j J'm 



'df' 


{dyjx=x. 


y=y 


m 


^9yVS 

j J'm 



i = 1, 2, . . . , n 
j = 1, 2, . . . , m 


(B8) 


(B9) 


(BIO) 


(BID 
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FORTRAN IV Subprograms for SPLINE Surface Fit 


Subroutine SPLIN2(X, Y, Z, N, M) : 

X Array of independent variable in increasing order 

Y Array of independent variable yj in increasing order 

Z Two-dimensional array of dependent variable f(Xj^,yj) corresponding to X and Y 
arrays 

This array must be dimensioned (20,20) in the calling program. 

N Number of values in X array; N < 20 

M Number of values in Y array; M < 20 
Purpose: 

To set up arrays for calling subroutine SPLINl and to store results from SPLINl for 
later use in the interpolation function subprograms. 

Program called: 

SPLINl 

Labeled COMMON: 

/BNDRY/BX1(20) , BXN(20) , JX, BY1(20) , BYM(20) , JY 

This labeled COMMON allows the user to specify one of the four sets of boundary 
conditions given by equations (B8) to (Bll): 



BX1(J) = < 


dxjx=x^ 

y=y 


j 




JX = 1 

J = 1, M 

j = 1, 2, . . . , m 

JX = 2 



JX = 1 

J = 1, M 
j = 1, 2, . . 

JX = 2 


m 
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BYia) = < 


1 — I 
by/x=x. 

y=Yi 


^ y=Yi 


r 


Vay/x=x. 

y=y. 


m 


BYM(1) 




'^sy 


2/x=x. 

y=y 


JY= 1 


JY = 2 


JY = 1 


JY = 2 


m 


I = 1, N 

i = l, 2, . . 


I = 1, N 

i = 1, 2, . . 


n 


n 


If the user does not specify any boundary conditions, the subroutines assume that the 
second particles are zero; that is, the arrays in labeled COMMON/BNDRY/ are ini- 
tialized to zero, and JX and JY are initialized to 2 by the block data subprogram. 
/SPUN/N1,M1,X1(20), Yl(20), Zl(20,20), ZX(20,20), ZY(20,20), ZYX(20,20) 

This labeled COMMON is used to transmit data from subroutine SPLIN2 to the func- 
tion subprograms. 

FORTRAN listings: 

$IBFTC SPLIN2 

SUBROUTINE SPL IN2 ( X , Y , Zt N, M ) 

DIMENSION X(?0),Y(20),Z(20,20>,ZI(20),ZP(20> 

COMMON /SPUN/ Nl, Ml, XH20), Yl (20), ZU 20, 20) fZX( 20, 20), 7 YI 20, 20 )» 

* ZYX(20,20) 

COMMON /8N0RY/ 8X 1( 20 ) ,BXN I 20 ) , JX ,8 Yl { 20 ) , OYM ( 20 ) , J Y 
N1 = N 
Ml = M 

DO 10 1=1, Nl 

10 XII I) = x( I ) 

DO 11 J=l,Ml 

11 YKJ) = Y(J) 

DO 12 1=1, Nl 
on 12 J=l,Ml 

12 Zl( I,J3 = Z(I,J) 

DO 3 I = l,Nl 

DO 1 J=1,M1 

1 ZI ( J) = Z 11 I, J ) 

CALL SPL INK Y1,ZI ,M1,ZP,RY1 (I ) ,8YM( I ) , JY) 

DO 2 J=l,Ml 

2 ZYI I, J ) = ZP( J ) 
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3 CONTINUE 

no 6 J = lfMl 
DO 4 I = ltNl 

4 ZI( I) = Zlf I,J) 

CALL SPL1NUX1#ZI,NI,ZP,BX1(J) tBXNi J) , JXJ 
DO 5 I=l,NI 

5 ZX{ I,J) = ZP(I) 

6 CONTINUE 

on 9 j=1tM1 
DO 7 1 = 1, N1 

7 Zim = ZYUfJ) 

CALL SPL IN 11X1,71 ,NI,ZP,0. 0,0*0,?) 

00 8 1=1, N1 

8 ZYX( I,J) = ZP( I ) 

9 CONTINUE 
RETURN 
END 

SIBFTC BOATA 

BLOCK DATA 

COMMON /BNDRY/ BX H 20 ) , B XN ( 20 ) , JX , B Y1 { 20) , BYM I 20 ) , J Y 
DATA BXl , BXN ,BY1 ,BYM, JX, JY/80#0. 0,2^2/ 

END 

Subroutine SPLINl: 

Purpose: 

This subroutine is called by SPIJN2 and is used to fit single spline curves through 
the data points on the grid lines. 

Calling program : 

SPUN2 

Labeled COMMON: 

None 

FORTRAN listing; 

$I8FTC SPL INI 

SUBROUTINE SPL IN 1 (.X, Y , N ,0 Y , B1 , BN, J ) 

DIMENSION X( I ) ,Y( I ) ,DY( I ) ,A (20) ,B(20) ,C( 20 ) ,01 20) 

Nl = N-1 
GO TO (4,5),J 

4 B( I ) = 1,0 
C( I ) = 0.0 
0( 1) = 61 
A(N) = 0..0 
8(N) = 1.0 
0(N) = BN 
GO TO 6 

5 BC I) = 2.0 
C( 1) = 1.0 

0(1) = 3.0«( Y( 2)-Y( 1) ) /( X( 2 )-X( 1) )-0.5«( X(2)-X( 1 )] *81 
A(N) = 1.0 
8(N) = 2.0 

DIN) = 3.0*( Y(N)“Y(N1) )/(X(N)-X(Nl) )+0.5*(XCN)-X(Nl))*BN 
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6 00 1 I=2fNl 

A( U = X( I+U-X( I ) 

B( I ) = 2.0*{ X( I +1 )-X( I-l) ) 

C( I ) - X( I )-X( I-l) 

1 0(1) = 3.0*( Y( I + l >*C( I )<=*24-Y(I )*( A( 1) **2-C( I )*<'2)-YU-l)*A( I )**2)/ 
♦ (C(I)*A(D) 

D( 1) = D( 1 )/fi( 1) 

DQ 2 1=2. N 

8(1) = 8( I )-A( I)«C(I-l)/B( I-l) 

2 0(1) = (D< I)-A( I )*0( I-l) )/8(I ) 

DY(N) = 0(N) 

DO 3 I=l,Nl 
K = N-I 

3 OY(K) = D( K)-C (K )*0Y (K+ 1) /B (K) 

RETURN 

END 


FUNCTION F(X1, Yl), FUNCTION FX(X1, Yl), FUNCTION FY(X1, Yl), and 
FUNCTION FXY(X1, Yl): 


XI The independent variable x. X(l) < XI < X(N) 

Yl The independent variable y. Y(l) < Yl < Y(M) 

Purpose: 

These four functions give interpolated values for the surface f(x, y) and the following 

p 

partial derivatives 0f/0x, 0f/0y, and 8 f/dx8y, respectively. 

Labeled COMMON: 

/SPUN/ 

FORTRAN listings : 

SIBFTC F 

FUNCTION F(X1,Y1) 

COMMON /SPIIN/ N.M, X( 20) ,Y( 20) ,Z( 20,20) , ZX(20. 20) ,ZY( 20,20) , 

* ZYX(20,20) 

G( F I, F II, OF I, OF 11, DEL, S) = (2. 0^= ( F I -F 1 1 ) +DFL* ( OF I + DF 1 1 ) ) ♦ ( S/ OEL ) <'* 

* 3 + ( 3.0*(F 1 1-F I )-0EL*( 2.0*DFI+DFI I ) ) #( S/OEL) **2 + 0FI*S + FI 
IF(Xl .LT.X( I ) ) STOP 

DO I 1=2, N 
K = I-l 

IF( XI .LE.X( I ) ) GO TO 2 

1 CONTINUE 
STOP 

2 IF( Yl ,LT.Y( 1 ) ) STOP 
DO 3 J=2,M 

L = J-1 

IF(Y1.LE.Y(J )) GO TO 4 

3 CONTINUE 
STOP 
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4 DEL = X(K+1)-X(K> 

S = Xl-X(K) 

FXJ = G(ZIK,L),Z(K+1,L) ,ZX(K,L) ,ZX(K+1,L),0EL,S> 

FXJl = G( Z(KtL+l) tZ{K+l,L*l) ,ZX1K,L+1) ,ZX(K+l,L+l) ,OEL,S) 

FYXJ = G( ZY(K,U,ZY(K+/1»L) ,ZYX(K,L) ,ZYX(K+l,L» ,DEL,S) 

FYXJl = G( ZY(K,L+l),ZYlK+ltL+l) »ZYX(K,L+U ,ZYX(K + l,L+n ,DEL,S) 

DEL = Y(L + U-Y(L» 

S = Yl-YIL) 

F = GJ FXJ, FXJl, FYXJ, FYXJl fOELfSJ 

RETURN 

END 

$IBFTC FX 

FUNCTION FX(X1,Y1) 

COMMON /SPLIN/ N,M,X120) ,Y(20) ,Z(20 ,20) ,ZX(20,20) ,ZY(?0,20), 

* ZYX( 20,20) 

G( Fr,FIl,OFI,OFIl,DEL,S) = ( 2 . 0* ( F I -F I I ) +OEL* ( OF H-OF 1 1 ) ) ♦ ( S /DEL ) ♦* 

* 3 + ( 3.0*(FI1-FI )-DEL*(2.0*DFH-0FI 1))*{S/0EL)«*2 + DFID'S + FI 
DG( FI,FII,0FI,DFI1,DEL,S) = ( 6. 0* ( F I-F 1 1 ) + 3 . 0*0EL* ( OF I+DF I 1 ) ) * 

* I S**2/0EL** 3 ) + l 6.0*( F I l-F I )-2. 0*DEL* ( 2. 0*DFI +UF I I ) )*( S/DEL«'*2 ) * 

* DFI 

IF( XI ,LT .X( 1 ) ) STOP 
DO 1 1=2, N 
K = I-l 

IFfXl,LE.X( I ) ) GO TO 2 

1 CONTINUE 
STOP 

2 IF( Y1 .LT.Y ( 1 ) ) STOP 
DO 3 J = 2,M 

L = J-1 

IF ( Y1 .LE.Y( J ) ) GO TO 4 

3 CONTINUE 
STOP 

4 DEL = X(K+1)-X{K) 

S = X l-X (K ) 

FXJ =0G( Z ( K,L ) , Z( K+UL ) ,ZX (K,LI , ZX( K+l ,L) ,OEL ,S ) 

FXJl =DG(Z(K,L+l),Z(K + l,L+l),ZX(K,L + l),ZX(K<-l,L4-l),0EL,S) 

FYXJ =DG(ZY(K,L),ZY(K+l,L), ZYX(K,L ) ,7YX(K+ l,L ) ,OEL ,S ) 

FYXJl =DG( ZY(K,L<-l),ZY(K+l,L+l),ZYX(K,L+l),ZYX(K + l,L+l),DEL,S) 

DEL = Y(L + 1)-Y(L) 

S = Yl-Y(L) 

FX =G1 FXJ, FXJl, FYXJ, FYXJl, DEL, S) 

RETURN 

END 

$I8FTC FY 

FUNCTION FY(XI,Y1) 

common /SPLIN/ N,M,X( 20) ,Y( 20) ,Z( 20 ,20) , ZX(20,20) ,Z Y( 20,20 ) , 

* ZYX( 20,20) 

G( FI,F Il,l)FI,DF 1 1 ,DFL , S) = 12. 0* ( FI -F 1 1 ) +DEL* ( OF I+DF I I ) )*(S/DEL ) ** 

* 3 + ( O.O’i'IFIl-FI )-0FL*(2.0*DFH-DFI 1))*(S/DEL)**2 * DFI*S + FI 
DGl FI,Fll, DFI,t)FI l,DFL,S) = I ft. 0* I F I-F 1 1 ) + 3. 0*DEL* ( DFI+OFIl) )# 

* ( S**2/DEL **3) + ( 6.0*( FI l-FI )- 2 . 0*DEL*( 2. 0*0FI +DF 1 1 ) )* ( S/DEL**2 ) + 

* DFI 

IFIXl .LT.Xm ) STOP 
00 1 I=2,N 
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K = 1-1 

lF(Xl.LF.Xi I ) ) GO TO 2 

1 CONTINUF 
STOP 

2 IF( Y1 .LT.YC H ) STOP 
00 3 J=2»M 

L = J-l 

IFiYl .LF.Y( J ) ) GO TO 4 

3 CONTINUE 
STOP 

4 DFL =: X(K+l)-X(K) 

S — X 1-“X ( K I 

FXJ = G(Z(K,L),Z(K+l,L»tZX(K,L> ,ZXIK+1,U,0EL,S) 

FXJl = G( Z (K ,L-H ) t Zl K + l»L+ 1 ) » ZX (K,L + 1 ) f ZX( K+l ,L + l ) »DEL»S) 

FYXJ = G(ZY(KtL),ZY(K+l,L) ,7YX(K,L)fZYX(K«^l,L) tOELfS) 

FYXJl = G(ZY{K,L + 1),ZY(K + 1,L+1) ,ZYX(K,L+1) »ZYXIK4-1,L+U ,D£L,S) 
DEL = Y(L+l)-Y(L» 

S = Yl-Y(L) 

FY =DG( FXJ ,FXJ 1»FYXJ,F YXJl ,OEL» S) 

RETURN 

END 

SIBFTC FXY 

FUNCTION FXY(X1,Y1) 

COMMON /SPLIN/ Nf MtX( ?0I f Y I 201 , Z( 20 t20) f ZX(20 ,20) f Z Yl 20 ,20 ) » 

* ZYX( 20, 20) 

DGIFI ,FIl,OFI,DFIl,DEL,S) = ( 6. 0* ( F 1-F 1 1 ) + 3. 0 *OEL* ( OF H-DF 1 1 )) * 

* (S**2/0EL**3) + lt».0*(FIl-FI )-2. 0*0FL* ( 2. 0*0F I+OF ID )*(S/DEL**2M- 

* DFI 

IFIXl .LT.X( D ) STOP 
DO 1 I=2,N 
K = I-l 

IFIXl. LE.X( I) ) GO TO 2 

1 CONTINUE 
STOP 

2 IF( Y1 .LT.YI I I ) STOP 
on 3 J=2,M 

L = J-1 

IF I Yl .LF.YI J ) ) GO TO 4 

3 CONTINUE 
STOP 

4 DEL = X(K+1)-X(K) 

S = Xl-XIK) 

FXJ =0G(Z(K,L),ZIK<^1,L),ZX(K,L) ,ZX(K«^1,L),DEL,S) 

FXJl =DG(Z(K,L+l),Z(K+l,L+l),ZX(K,L+l),ZX(K+l,L+l),OEL,S) 

FYXJ =DGI ZY(K,L) ,7Y(K+l,L) ,ZYX( K,L) ,ZYX(K+1,L) ,OEL,S) 

FYXJl =DG(ZY(K,L+n ,ZY(K+l,L+D ,ZYX(K,L«^1) ,Z YX I K+l , L+U ,DEL , S > 
DEL = Y(L + 1)-Y(U 
S = Yl-YIL) 

FXY = DGl FXJ, FXJl, FYXJ, FYXJl , DEL, S) 

RFTURN 

END 
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Program use: 

The user must call subroutine SPLIN2 with the proper arguments only once for the 
surface to be fit. After this call the user may use any or all of the function sub- 
programs . The user must make sure that the arguments used in the functions are 
within the ranges of the original rectangular grid. These programs should not be 
used for extrapolation. 
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APPENDIX C 


CONVERGENCE EXPERIMENT 
Problem Description 


Given the equation 


y”(x)=f(x,y) (Cl) 

it is desired to investigate the nature of the convergence by Picard iteration to the solu- 
tion y*(x)when f(x,y) is, for a given x, a normally distributed random variable. In 
particular, it is desired to determine the effect of the standard deviation of f(x,y), at a 
given X, on the convergence. 

A general analytical analysis has not been attempted. Instead a particular expression 
has been chosen for f(x,y). Solutions of equation (Cl) have been obtained numerically by 
the Clenshaw-Norton method (ref. 2) for different standard deviations, and the results 
compared. 


Numerical Experiment 

Let f(x,y) = 4y + 6(x), where 6(x) is, for a given x, a normally distributed ran- 
dom variable with zero mean and standard deviation <^. Hence, f(x,y) is, for a given x, 
a normally distributed random variable with mean 4y and standard deviation a^. The 
eouation 

y”(x) = 4y^ + 5(x) 0 < x < 1 (C2) 

with the boundary conditions y(0) = 0 and y(l) = 1 was numerically solved for the fol- 
lowing values of the standard deviation Og: 0, 0.05, 0.1, 0.2, 0.5. Note that the 
Clenshaw-Norton method evaluates the right side of equation (C2) only at the N argu- 
ments of an N - 1 degree Chebyshev curve fit 

9 N-1 

4y^ + 6(x) = 

k=0 


where Tj^(x) are the Chebyshev ploynomials of degree k. For this experiment, N = 17. 
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Results 

The relative rate of convergence is shown in figure 37 for the various standard de- 
viations. The curve labeled ctq = 0 is, of course, the exact solution. The error in fig- 
ure 37 is defined as the maximum difference in two successive iterations of the corre- 
sponding coefficients a^^ (eq. (C3)). For Og > 0, a general tendency to oscillate about a 
mean error is noted as the number of iterations increase; the magnitude of this mean 
error is proportional to Og, as might be expected. 

The results of figure 37, while interesting, are not of major significance in relation 
to the problem encountered in this report. Of greater interest is the relative error of 
the solution. Since equation (C2) is solved as a boundary value problem, the relative 
error was investigated at x = 0. 5. The results are given in table IV. 
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TABLE IV. - EFFECTS OF STANDARD DEVIATION 


(a) On relative error during convergence 


Iteration 

Standard deviation, Og 


0 

0.5 

0.2 

0.1 

0.05 


Relative error at 

X = 0. 5, 

(y(x) - y» 

(x))/y*(x) 

1 

-0.685 

-0.685 

-0.685 

-0.685 

-0.685 

2 

.235 

.189 

.217 

.226 

.231 

3 

-.098 

-.098 

-.098 

-.098 

-.098 

4 

,033 

-.042 

.003 

.018 

.026 

5 

-.012 

.154 

.054 

.021 

.005 

6 

.004 

-.094 

-.034 

-.015 

-.005 

7 

-.001 

.012 

.004 

.001 

-.000 

8 

.000 

-.034 

-.013 

-.006 

-.003 

9 

-.000 

-.093 

-.037 

-.019 

-.009 

10 

.000 

-.052 

-.021 

-.010 

-.005 

11 

-.000 

-.042 

-.017 

-.008 

-.004 


(b) On standard deviation a ^ 


Iterations 
used in 
computing 

Standard deviation, 

0 

0.5 

0.2 

0.1 

0.05 


Standard deviation of y(x) at x = 0. 5, 

2 to 11 
4 to 11 

0.010 

.002 

0.013 

.011 

0.010 

.004 

0.010 

.002 

0.010 

.002 


Discussion of Results 

The effect of increasing the standard deviation Og is shown in table IV. The rela- 
tive error tabulated in table IV(a) is taken with respect to the ’’true” value of y(0. 5) 

(i.e., with ff(x) = 0 in eq. (C2)). Of greatest interest is the standard deviation, of 
y(0. 5) about the true value, as shown in table IV(b). There, is given, calculated 
from two choices of sample iterations; iterations 2 to 11 and iterations 4 to 11. It is ob- 
served that when the last 10 iterations are employed in the calculation of Oy, is in- 
sensitive to CTg. This observation tends to corroborate the method employed in this re- 
port of averaging successive iterations. 

Averaging also accomplishes another economy. The evaluations of f(x,y), with the 
associated random errors, at the Chebyshev arguments x. can be likened to experi- 
mentally obtained data at this same argument. In the case of experimental data, the usual 
procedure would be to use a least-squares fit. This could be accomplished in the present 

110 


I M Ml I I I ■■!■■■ ■ I 




“ I I 


I I 


I 



program by simply truncating the expansion of equation (C3) at some M < N. A test was 
made by setting M = 5. While the resulting fit was smoother (after each iteration) the 
results, as presented in table IV, were not affected. 


Conclusions 

This numerical experiment on the convergence of equation (C2) is not necessarily ex- 
pected to apply to all second-order differential stochastic equations. In fact, there is no 
assurance that this experiment accurately represents the problem in ENEC where 

f(x,y) = C • n(x) (C4) 

This experiment was primarily intended to give some insight into the problem. With this 
in mind it must be noted that, by the Central Limit Theorem (appendix A of ref. 1), 
n(x) + 6(x) approaches n(x) as the number of electron histories becomes very large. 
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APPENDIX D 


I/VIPROVED SQUARE ROOT ROUTINE 

standard computer library subroutines for the calculation of square roots have accu- 
racy greater than is needed for Monte Carlo work. Sacrificing some of this accuracy for 
a faster square root subroutine has proved to be helpful in increasing the efficiency of 
Monte Carlo codes. 

Three modifications were made to the IBM FORTRAN IV version 13 library square 
root subroutine to effect the increase in speed: 

(1) The number of Newton-Raphson iterations, applied to the starting approximation, 
was reduced from three to two or one. 

(2) The indexing on the iteration loop is removed, and the iteration equation was 
written twice or once. 

(3) The starting approximation was changed. 

The last modification was necessary to maintain good accuracy with the reduced num- 
ber of iterations. The starting approximation used is given as equation (11”) in refer- 
ence 11. The modifications added several words to the subroutine, but it was felt that the 
increase in speed more than compensated for this. 

Table V lists results from tests run on square root subroutines with various com- 
binations of the three modifications. Each test consisted of 100 000 arguments exponen- 

38 

tially distributed over the range from 0 to 10 and was conducted on an IBM 7094 com- 
puter. 

Test 1 is the library square root. Test 2 is the starting approximation used in the 
library square root with two iterations and no indexing. Test 3 is the same as test 2 but 

TABLE V. - SQUARE ROOT SUBROUTINE TEST RESULTS 


[F and G represent true and test square roots, respectively.] 


Test 

Ratio of library 
to test square 
root time 

/ ^ 
In/ 

\l/2 

F-Gin 

^7 

MAxl^ " 

1 F 1 

1 

1.0 

3. IxlO'® 

7.25x10"® 

2 

1.65 

3.65x10"'^ 

1.5x10"® 

3 

2.41 

5. 8X10"^ 

1.7x10'® 

4 

1.77 

9.7X10"® 

2.0x10"'^ 

5 

2.55 

3.7X10"'* 

6.4x10"'* 
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with one iteration. Tests 4 and 5 use the more accurate starting approximation of the 
reference with two and one iterations, respectively, without indexing. 

As a final test, the square root subroutine that gave a time ratio of 1.77 was tried in 
ENEC which used the square root subroutine 44 percent of the running time. The saving 
in computer time for this code was considerable, while the results compared favorably 
with those using the slower library subroutine. The listing of FSQR, the subroutine used 
in this final test, is given in appendix E. 
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APPENDIX E 


MACHINE LANGUAGE SUBROUTINES 
RANDOM 


Purpose ; 

A double -entry subroutine for generating pseudorandom numbers within the range 
from 0 to 1. 

Calling sequence: 

CALL SAND(RO) 

CALL RAND(R) 


Use: 

The first call to this subroutine must be CALL SAND(RO). This call is necessary to 
set up addresses in RAND, The statement CALL RAND(R) will cause the next 
pseudorandom number to be generated. The normalized floating point number is 
stored in R, while the fixed point number is stored in RO. 

Method: 

The pseudorandom number sequence is generated by the low order 36 bits of the 
product where K = 5^^, r^_j^ is the previous pseudorandom number, and 

r^ = 1 (see ref. 4). This fixed point number is then floated and normalized so that 
the result is in the range from 0 to 1. 

Remarks : 

This method is dependent on the computer word length and was specifically designed 
for the IBM 7094. 

Map listing: 


$ IBMAD 

random 






ENTRY 

RAND 





ENTRY 

SAND 




SAND 

CLA 

3, A 





sta 

B 

MULTIPLIER 

IN RANDOM NO. 

GENERATOR 


CLA 

ONE 

SET DAM = 

TO 1 FOR FIRST 

RANDOM NO 

DAW 

STO* 

3»A 





TRA 





rand 

SAVE 

(4) 





LDO* 

B 





MPY 

CONS 

BY DAM 



B 

STO 


STORE THE 

LOW ORDER PART 

AT DAM 


CLA 

flc 

FLOAT normalize * AND 



LLS 

27 

ROUND the 




FAD 

C 

random no. 



R 

STO* 

3*4 





TRA 

1 *4 




FLC 

OCT 

000000U00200 

EXPONENT OF RANDOM NO. 
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CONS 

DEC 

30517578125 

5 EXP 15 

ONE 

DEC 

1 

ONE 

C 

OCT 

END 

170000000200 

NORMALIZING CONSTANT 


BCREAD 

Purpose: 

Subroutine BCREAD allows the programer to read absolute binary cards with a maxi- 
mum of 22 words per card. BCREAD is to be used in conjunction with BCDUMP. 
Calling sequence: 

CALL BCREAD(A,B), where A is the first data word to be read, and B is the last 
data word to be read. If the address of A equals the address of B, one word is 
read. 

Remarks: 

The address of A must be less than or equal to the address of B. BCREAD makes 
use of the file definition subroutine . READS. 

Map listings: 


SIBMAP 

BCREAD 




TTL 

BCREAD SUBROUTINE FOR IBSYS 


lbl 

BCREAD 



ENTRY 

BCREAD 


BCREAD 

SAVE 

l*2t4 



CLA 

3,4 

PICK UP THE FIRST ARGUMENT 


LDO 

4,4 

PICK UP THE SECOND ARGUMENT 


TLQ 

*-1-2 

MAKE SURE THE LARGEST 


XCA 


ARGUMENT IS IN THE AC 


STO 

TEMP 



SUB 

TEMP 



PAX 

0,1 

PUT WORD COUNT + 1 


TXT 

*+l ,1 ,1 

INTO INDEX 1 


LXA 

TEMP, 2 

PICK UP THE FIRST LOAD ADDRESS 


SXA 

1X1,1 



SXA 

1X2,2 



CLA* 

INS 



STA 




TSX 

•CLOSE, 4 



MON 




CLA* 

READS 



STA 

MON 



STA 

READ2 



STA 

SHUT 



TSX 

•OPEN, 4 


MON 

MON 

** 


IXl 

AXT 

**,1 

HOLDS THE WORD COUNT 

1X2 

AXT 

**,2 

HOLDS THE LOADING ADDRESS 

SXA 

SXA 

10,2 



txl 

LASTCfl ,22 


READ 

TSX 

•READ, 4 


READ2 

PZE 

**, ,E0B 
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10 

PZE 

lOCPN 

lOCD 

LASTC 

TXI 

TXI 

TRA 

CLA 

DONE 

STO 

SXD 

TRA 

TRA 

SHUT 

AXT 

SXA 

AXT 

SXD 

TSX 

MON 

EOB 

RETURN 

CALL 

ERR 

TRA 

call 

EOF 

CALL 

ERR2 

PZE 

E0B2 

PZE 

INS 

PZE 

READS 

PZE 

TEMP 

PZE 


END 


EOF* »ERR 
**» »2 
**» *22 
*+l»l,-22 
*+l»2*22 
SXA 
DONE 
*-2 
10,1 
READ 
*+l 
SXA,4 
LASTC-1 ,4 

22.4 

10.4 

.CLOSE, 4 
** 

BCREAO 

.FXEM.(E0B2) 

EOF 

.FXEM. (ERR2) 
EXIT 

35 

36 

.UN05, 

.READ5 


SKIP first word and CHECKSUM 


REDUCE THE WORD COUNT 


SIBMAP .READS 
ENTRY 

.READS PZE 
READS file 
END 


• READS 
READS 

,IN1,READY,INPUT,BLK=28,MULTIREEL,MXBIN,N0LIST 


BCDUMP 


Purpose: 

Subroutine BCDUMP allows the programer to punch out data in an absolute binary 
format with a maximum of 22 words per card. BCDUMP is meant to be used in con- 
junction with BCREAD. 

Calling sequence : 

CALL BCDUMP (A, B,K), where A is the first data word to be punched and B is 
the last. If the address of A equals the address of B, one word is punched. 

K controls card numbering. If K equals zero or is missing, each call to BCDUMP 
will start numbering cards with 000. If K is not equal to zero, the numbering con- 
tinues in sequential order starting with 000. 

Remarks: 

The address of A must be less than or equal to the address of B. BCDUMP makes 
use of the file definition subroutine . PCH. . 
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Map listings: 


SIBMAP 

BCDUMP 




ttl 

BCDUMP ROUTINE 

FOR IBSYS 


ENTRY 

BCDUMP 


BCDUMP 

SAVE 

1»2*4 



CLA 


IS THERE A 


PDX 

0,2 

THIRD 


TXL 

*+2,2,2 

ARGUMENT 


N2T* 

5,4 

YES. IS IT = 0 


SXA 

CNUMtO 

YES 


CLA 

3,4 

PICK UP FIRST ARGUMENT 


LDO 

4,4 

PICK UP SECOND ARGUMENT 


TLO 

*+2 



XCA 




STQ 

WDl 

WOl HAS THE FIRST ADDRESS 


LXA 

WD1,1 

FIRST LOCATION IN INDEX 1 


SU9 

WDl 



PAX 

0,2 

THE NO. OF WORDS OUTPUTED IN INDEX 2 


TXI 

*+1,2,1 

TRUE WORD COUNT 


SXA 

1X1,1 



SXA 

1X2,2 



CLA* 

OUT 



STA 

RITE+1 



STA 

MON 



STA 

CLSE 



PAX 

0,1 



TXI 

*+1,1,1 



SXA 

*+1,1 



LDI 




LNT 

040000 



TRA 

*+2 



TRA 




TSX 

.0PEN,4 


MON 

MON 



IXl 

AXT 

**, 1 


1X2 

AXT 

**, 2 


TEST 

TXL 

LASTC,2,22 

ONLY 22 WORDS OR LESS LEFT 


TlX 

*+1,2,22 



SXA 

1X2,2 



AXT 

22,2 


TESTA 

TXI 

*+1,2,320 



SXO 

WDl ,2 



TIX 

*+1,2,320 



SXA 

CLA,1 



SXA 

WDl ,1 


LOOP 

TXI 

*+1,1,22 



SXA 

1X1,1 



AXT 

23,4 


CLEAR 

ST2 

CKSUM+23,4 

CLEAR THE BUFFER 


Tlx 

*-1,4,1 



AXT 

0,4 


CLA 

CLA 

**,4 

FILL THE BUFFER WITH 


STO 

CKSUM+1,4 

NEW DATA 


TXI 

*+l ,4,-1 



TIX 

*-3,2,1 


CNUM 

AXT 

**,1 

CONSECUTIVELY 


CLA 

HUNBIT 

NUMBER 


ARS 

1 

THE ' / 
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I 


I 




TXL 

♦+2flf99 

BCDUMP 


TXI 

*-2tlf-100 

CARDS 


STA 

GP+1 

from 


CLA 

BITT 

ZERO 


ARS 

1 

TO 


txl 

♦•f2»l*9 

999 


TXI 

*-2*l*-10 



STO 

W0RD3 



CLA 

BITU 



ARS 

1 



TXL 

*+2f IfO 



TXI 

*—2 » 1 f — 1 



ORS 

W0RD3 



LXA 

CNUM*1 



TXI 




txl 

*+2*1 >999 



AXT 

c*i 



SXA 

CNUM*1 



AXT 

22*1 

COMPUTE 


cal 

WDl 

THE 


ACL 

GP*1 

CHECK SUM 


Tlx 

*•-1*1*1 



SLW 

C<SUM 


RITE 

TSX 

•WRITE*4 

WRITE THE BINARY CARD ON THE OUTPUT TAPE 


PZE 

*♦* *EOF 



IDCT 

WD1**28 

CHANGE TO TRANSMIT FOR DIRECT COUPLE LMLR 


TRA 

IXl 


RETURN 

TRA 

*+l 



AXT 

1X1*1 



SXA 

RETURN-1*! 



TSX 

.CLOSE*4 


CLSE 

MON 




RETURN 

BCDUMP 


lastc 

CLA 

RETURN 



STO 

RETURN-1 



TRA 

TESTA 


EOF 

TSX 

.FXEM.*4 

ERROR ON END OF FILE LMLR 


TXI 

*+3* *1 

LMLR 


PZE 

** ,OUT+l 

LMLR 


PTH 

E0F2**2 

LMLR 


call 

EXIT 


E0F2 

BCI 

2*EOF BCDUMP 


WDl 

PZE 



CKSUM 

BSS 

23 


GO 

OCT 

420041004040 



OCT 

104020400000 


WORD3 

PZE 




PZE 



HUNBIT 

OCT 

2000 


BITU 

OCT 

20000000 


BITT 

OCT 

200000000000 


OUT 

PZE 

• PCH. 



END 



$IBMAP 

.PCH. 




ENTRY 

• PCH. 


• PCH* 

PZE 

PCH 


PCH 

file 

*PP * READY* OUTPUT *BLK=28 *MULT I REEL *B I N * NDL 1ST 


END 
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FSQR 


Purpose: 

A double -entry subroutine to compute the square root by the method described in 
appendix D. 

Calling sequence: 

SQRT(X), ASQRT(X) 

Remarks: 

The entry point SQRT computes the square root of the argument and returns with the 
result in the accumulator. The entry point ASQRT forces the sign of the argument 
positive before computing the square root. 

Map listing: 


SIBMAP 

FSOR 



ENTRY 

ASQRT 


ENTRY 

SORT 

ASQRT 

CLA* 

3*4 SQUARE ROOT SUBROUTINE 


T2E 



SSP 

EVE-S FIRST APPROXIMATION (11 


tra 

BEGIN TWO NEWTDN-RAPHSON ITERATIONS 

SORT 

CLA* 

3*4 03/08/66 


T2E 

1*4 


TMI 

ERROR 

BEGIN 

STO 

BUFF 


ana 

K1 


TZE 



SUB 

K3 


ADD 

BUFF 


AR5 

1 


ADD 

K2 


STO 

BUFF+1 


CLA 

BUFF 


FDP 

XCA 

BUFF+1 


FAD 

BUFF+1 


SUB 

K1 


STO 

BUFF+1 


CLA 

BUFF 


FDP 

XCA 

BUFF+1 


FAD 

BUFF+1 


SUB 

Kl 


TRA 

1*4 

ERROR 

SXA 

SYSL0C*4 


SXA 

LINK*4 


CALL 

•FXEM. (ESORT) 


ORG 

*^1 

ESORT 

MTW 

• SORTN* *3 


LXA 

SSP 

LINK*4 


TRA 

BEGIN 
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K1 

OCT 

001000000000 

K2 

OCT 

100356300000 

K3 

OCT 

000356300000 

BUFF 

BSS 

2 

L INK 

LDIR 

END 



APPENDIX F 


SAMPLE PROBLEM 

The example given in this appendix is ENEC's solution to the equation 

(p”{x) = 50n(x) 

for <p’(0) = -13 and ^(0) = 0, where <p(x) is the potential distribution and n(x) is the 
dimensionless electron density distribution. The scattering gas is to be argon at a tem- 
perature of 1000^ K. The reduced pressure in torr times the interelectrode spacing in 
centimeters is taken to be 2. The number of particles for each iteration is 500, 

and averaging is to take place over 10 iterations after convergence to 0. 8 is achieved. 

The Chebyshev fit is to use 17 data points and the initial fit of ^(x) was obtained from a 
previous ENEC run. 

Tables VI to XII give the input and sample printed output of the programs that con- 
struct the tables needed by ENEC (see section Preparation of Input Tables). The input 
data were obtained for argon. To get a three dimensional perspective for the argon elec- 
tron differential scattering cross-section surface u(0,E), the surface was drawn on a 
digital mechanical plotter and is depicted in figure 38. Output from ENEC (described in 
the section ENEC Output) for this example is given in figure 39 and table XTTT . This sam- 
ple problem was taken from reference 12. This reference contains many results obtained 
from the ENEC code. 
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TABLE VI. - FIRST PAGE OF PRINTED OUTPUT FROM PROGRAM CVEL -In R 


I VEL 


I 

7.6246 

6 . 5260 

6.01 52 

5.6787 

5.4274 

5.2267 

5.0597 

4.9166 

9 

4.7914 

4.6802 

4.5801 

4.4891 

4.4057 

4.3288 

4.2573 

4.1906 

17 

4.1281 

4.0693 

4.0137 

3.9611 

3.9110 

3.8634 

3.8180 

3.7745 

23 

3.7328 

3.6928 

3.6543 

3.5173 

3.5816 

3.5471 

3.5137 

3.4815 

33 

3.4502 

3.4199 

3.3905 

3.3619 

3.3342 

3.3071 

3.2808 

3.2552 

41 

3.2302 

3.2058 

3.1820 

3.1587 

3.1360 

3.1138 

3.0920 

3.0707 

49 

3.0499 

3.0295 

3.0095 

2.9899 

2.9707 

2.9518 

2.9333 

2.9151 

57 

2.8972 

2.8797 

2.8624 

2.8455 

2.8288 

2.8124 

2.7963 

2.7604 

65 

2.7648 

2.7494 

2.7343 

2.7193 

2.7046 

2.6901 

2.6759 

2.6618 

73 

2.6479 

2.6342 

2.6207 

2.6073 

2.5942 

2.5812 

2.5684 

2.5557 

81 

2.5432 

2.5309 

2.5187 

2.5066 

2.4947 

2.4830 

2,4713 

2.4598 

89 

2.4485 

2.4372 

2.4261 

2.4151 

2.4043 

2.3935 

2. 3829 

2.3723 

97 

2.3619 

2. 3516 

2.3414 

2.3313 

2.3213 

2.3114 

2,3016 

2.2919 

105 

2.2823 

2.2728 

2.2633 

2.2540 

2.2447 

2.2355 

2.2265 

2,2174 

113 

2.2085 

2. 1997 

2 . 1909 

2.1822 

2.1736 

2.1650 

2.1566 

2.1482 

121 

2.1398 

2.1316 

2.1234 

2 . 115 ? 

2.1072 

2 . 099 ? 

2.0912 

2.0834 

129 

2.0755 

2.0678 

2.0601 

2.0525 

2.0449 

2,0374 

2.0299 

2,0225 

137 

2.0151 

2. 0078 

2.0006 

1.9934 

1.9863 

1.9792 

1.9721 

1.9551 

145 

1 .9582 

1.9513 

1.9444 

1.9376 

1.9309 

1.9242 

1.9175 

1.9109 

153 

1.9043 

1.8978 

1.8913 

1.8843 

1.8784 

1.8720 

1.8657 

1.8594 

161 

1.8532 

1.8470 

1.8408 

1.8347 

1.8286 

1.8225 

1.8165 

1.8105 

169 

1.8045 

1.7986 

1.7927 

1.7869 

1.781 1 

1.7753 

1.7695 

1.7638 

177 

1.7582 

1.7525 

1.7469 

1.7413 

1.7357 

1.7302 

1.7247 

1.7193 

185 

1.7138 

1. 7084 

1.7030 

1,6977 

1-6924 

1.6871 

1.6818 

1.6766 

193 

1.6714 

1.6662 

1.6610 

1.6559 

1.6508 

1.6457 

1.6407 

1.6357 

201 

1.6307 

1.6257 

1.6207 

1.6158 

1.6109 

1.6060 

1.6012 

1.5963 

209 

1.5915 

1.5867 

1.5820 

1.5772 

1.5725 

1.5678 

1.5632 

1.5585 

217 

1.5539 

1.5493 

1.5447 

1.5401 

1.5356 

1.5310 

1.5265 

1.5221 

225 

1.5176 

1. 5132 

1.5087 

1.5043 

1.4999 

1.4956 

1.4912 

1.4869 

233 

1.4826 

1.4783 

1.4740 

1.4698 

1.4655 

1.4613 

1.4571 

1.4529 

241 

1.4488 

1.4446 

1.4405 

1.4364 

1.4323 

1.4282 

1.4241 

1.4201 

249 

1.4160 

1.4120 

1.4080 

1.4040 

1.4001 

1.3961 

1.3922 

1.3882 


TABLE VII. - OUTPUT OF INPUT 
TO PROGRAM MFP a(E^/^) 


SQKT ( SV ) 

SIGMA-S 

0 . 

36.40 

0 . 1000000 

25.00 

0.2000000 

13,25 

0.3000000 

6 . 70 

0.4000000 

3.30 

0.5000000 

1 . 70 

0.6000000 

1.40 

0.7000000 

1.60 

0. 8000000 

2.35 

0.9000000 

3.40 

1. 0000000 

4.80 

1.4142136 

12.00 

1.7888543 

20.20 

2.0000000 

25.40 

2.2360680 

32.20 

2.5884358 

45.00 

2.8284271 

54.50 

3.0000000 

63 . 00 

3.2093612 

72.20 

3.5355338 

80.00 

3.8729833 

84 . 16 


TABLE Vm. 

- PRINTED OUTPUT OF PROGRAM 

TABLE Vin. 

- Continued. 

PRINTED OUTPUT OF 


MFP X(E) 



PROGRAM MFP X(E) 


EV 

MFP 

56 

0.555 

0.53687751 




57 

0.565 

0. 52264239 

1 

0.005 

0.03545675 

58 

0.575 

0. 50876888 

2 

0.015 

0.04547802 

59 

0. 585 

0.49530750 

3 

0.025 

0. 05670404 

60 

0.595 

0.48229624 

4 

0.035 

0.06943773 

61 

0.605 

0.46976049 

5 

0.045 

0.08314588 

62 

0,615 

0.45771582 

6 

0.055 

0.09728152 

63 

0.625 

0.44616985 

7 

0.065 

0.11181470 

64 

0.635 

0.43512379 

8 

0.075 

0. 12682515 

65 

0.645 

0.42457313 

9 

0.085 

0. 14252448 

66 

0.655 

0.41449207 

10 

0.095 

0.15929276 

67 

0.665 

0.40483951 

11 

0.105 

0.17744290 

68 

0.675 

0.39557862 

12 

0. 115 

0.19708641 

69 

0.685 

0. 38667697 

13 

0.125 

0.21831543 

70 

0.695 

0. 37810599 

14 

0.135 

0.24120887 

71 

0,705 

0. 36984034 

15 

0.145 

0.26582742 

72 

0,715 

0.36185761 

16 

0.155 

0.29220729 

73 

0.725 

0.35413782 

17 

0. 165 

0. 32035189 

74 

0.735 

0.34666324 

16 

0.175 

0. 35019431 

75 

0.745 

0. 33941 796 

19 

0.185 

0.38156812 

76 

0. 755 

0. 33238783 

20 

0.195 

0.41421016 

77 

0.765 

0. 32556009 

21 

0.205 

0.44774819 

78 

0.775 

0. 31892335 

22 

0.215 

0.48169303 

79 

0.785 

0. 31246728 

23 

0.225 

0.51543956 

80 

0.795 

0.30618261 

24 

0.235 

0, 54827990 

81 

0.805 

0. 30006094 

25 

0,245 

0.5 7943085 

82 

0.8L5 

0.29409487 

26 

0.255 

0. 60808366 

83 

0.8 25 

0.28828214 

27 

0.265 

0.63362976 

84 

0.835 

0. 28762451 

28 

0.275 

0.65575616 

85 

0.845 

0. 27712286 

29 

0.285 

0.67432405 

86 

0.955 

0.27177718 

30 

0.295 

n. 68935928 

87 

0. 865 

0.26658668 

31 

0.305 

0. 70103466 

88 

0.875 

0. 26154993 

32 

0.315 

0. 709641 84 

89 

0.885 

0.25666498 

33 

0.325 

0. 71555825 

90 

0.895 

0. 25192939 

34 

0.335 

0. 71921379 

91 

0.905 

0.24734043 

35 

0.345 

0. 72106143 

92 

0,915 

n. 24289504 

36 

0.355 

0.72155368 

93 

0.925 

0. 23858999 

37 

0.365 

0. 721 11481 

94 

0,935 

0. 23442183 

38 

0.375 

0. 71989438 

95 

0.945 

0. 23038704 

39 

0.385 

0. 71779265 

96 

0. 955 

0, 22648200 

40 

0.395 

0. 71471675 

97 

0.965 

0.22270304 

41 

0.405 

0.71059322 

98 

0,975 

0.21904650 

42 

0.415 

0. 7053691 1 

99 

0.985 

0.21550869 

43 

0.425 

0.69901259 

100 

0.995 

0.21208598 

44 

0.435 

0.69151314 

lOl 

1.250 

0. 15176453 

45 

0.445 

0.68288124 

102 

1.750 

0.09863517 

46 

0.455 

0.67314735 

103 

2.250 

0,07346389 

47 

0.465 

0. 66236047 

104 

2.750 

0. 05875729 

48 

0.475 

0.6505 86 06 

105 

3.250 

0.04919928 

49 

0.485 

0.63790367 

106 

3.750 

0.04247487 

50 

0.495 

0.62440746 

107 

4.250 

0.03735838 

51 

0. 505 

0.61027010 

108 

4.750 

0,033201 96 

52 

0.515 

0.59571803 

109 

5.250 

0.02966338 

53 

0.525 

0. 58094945 

110 

5.750 

0.02667541 

54 

0.535 

0.566131 75 

111 

6.250 

0.02423191 

55 

0.545 

0. 55140359 

112 

6.750 

0.02226726 




113 

7. 250 

0.02066121 




114 

7.750 

0. 01922455 




115 

8.250 

0.01785111 


123 



TABLE Vni. - Concluded. PRINTED OUTPUT OF 


PROGRAM MFP X(E) 


116 

8.750 

0. 01658369 

117 

9.250 

0 . 01553 L 71 

118 

9.750 

0.01469191 

119 

10.250 

0.01404402 

120 

10.750 

0.01356194 

121 

11.250 

0.01320282 

122 

11.750 

0.01292968 

123 

12.250 

0.01271599 

124 

12.750 

0 . 012541 88 

125 

13.250 

0.01239486 

126 


0.01199999 

* 0l + EXIT 

IN MFP 



TABLE IX - INPUT TO PROGRAM ARGON FOR ct(0, E) 


SDATA 





0 . 

15 . 

28 . 

43 . 

59 . 

74.5 

90 . 

105.5 

121 . 

137 . 

152 . 5167.5 

180 . 

• 

01 

2 

.73 

2.53 

2.38 

2.21 

2.05 

1.91 

1.79 

1.70 

1.62 

1.56 

1.52 

1.49 

1.49 

• 

05 

2 

.36 

1.96 

1.67 

1.38 

1.13 

.933 

.780 

.662 

.574 

.509 

.468 

.447 

.441 


.1 

2 

.06 

1.55 

1.20 

.880 

.624 

.444 

.316 

.228 

.169 

.130 

.107 

.096 

.093 


.2 

1 

.68 

1.07 

.687 

.386 

.188 

.081 

.027 

.005 

.000 

.002 

.006 

.008 

.009 


.3 

1 

.45 

.788 

.417 

.165 

.039 

.001 

.007 

.028 

.052 

.072 

.085 

.091 

.093 


.4 

1 

.31 

.613 

.262 

.062 

.001 

.017 

.060 

.102 

.133 

.153 

.162 

.166 

.166 


.5 

1 

.23 

.502 

.170 

.018 

.009 

.062 

.123 

.168 

.192 

.199 

.198 

.195 

.194 


.6 

1 

. 2C 

.431 

.116 

.003 

.033 

.108 

.173 

.207 

.214 

.204 

.190 

.180 

.177 

1 

.1 


.00 

.030 

.050 

.100 

.163 

.190 

.199 

.172 

.131 

.090 

.063 

.032 

.000 

2 

.0 


.06 

.100 

.170 

.300 

.390 

.420 

.359 

.280 

.180 

.115 

.120 

.100 

.07 

2 

.4 


.15 

.200 

.240 

.380 

.490 

.500 

.420 

.320 

.200 

.140 

.160 

.150 

.12 

2 

.8 


.28 

.300 

.300 

.450 

.570 

.590 

.470 

.350 

.216 

.160 

.200 

.210 

.18 

3 

.2 


.42 

.425 

.370 

• 540 

.660 

.670 

.530 

.390 

.240 

.188 

.240 

• 282 

• 26 

4 

.0 


.84 

.700 

.540 

.580 

.820 

.820 

.630 

.450 

.280 

.270 

.336 

.460 

#44 

5 

.0 

1 

.68 

1.10 

.770 

.850 

1.00 

1.00 

.750 

.530 

.340 

.410 

#660 

.800 

• 80 

6 

.7 

4 

.85 

2.42 

1.36 

1.13 

1.22 

1.18 

.915 

.620 

.420 

.660 

1.37 

1.88 

2.14 

8 

.0 

7 

.00 

4.45 

2.06 

1.18 

1.19 

1.16 

.885 

.640 

.458 

.824 

1.80 

2.68 

3.17 

9 

.0 

8 

.25 

5.75 

2.84 

1.60 

1.33 

1.19 

.947 

.690 

.563 

1.06 

2.10 

3.03 

3.50 

1C 

.3 

9 

.60 

7.15 

4.42 

2.27 

1.27 

1.11 

1.00 

.775 

.610 

1.16 

2.38 

3.30 

3.70 

12 

.5 

i : 

3.0 

9.30 

4.93 

2.38 

1.31 

.960 

.830 

.690 

.690 

1.23 

2.58 

3.49 

3.90 
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TABLE X. - OUTPUT FROM PROGRAM ARGON 


i 


cos e 


cr(c os~^t, E)dt 
(t(E) 


SIGMA 


EV/COSI THETA ) 

-l.O 

-0.8 

-0.6 

-0.4 

-0.2 

-0.0 

0.2 

0.4 

0.6 

0.8 

l.O 

0.01 

-0. 

0.0816 

0. 1661 

0.2536 

0.3446 

0.4392 

0,5378 

0.6417 

0.7517 

0-8693 

1.0000 

0-05 

-0. 

0.0520 

0. 1098 

0.1745 

0-2471 

0.3290 

0.4219 

0,5289 

0.6538 

0.8033 

l.OOOO 

0.10 

>0. 

0.0235 

0.0535 

0.0918 

0.1406 

0-2029 

0,2824 

0.3862 

0. 5223 

0.7081 

1.0000 

0.20 

-0- 

0.0041 

0.0052 

0.0053 

0.0079 

0.0192 

0.0466 

0,1094 

0,2289 

0.4601 

1.0000 

0.30 

-0. 

0.0870 

0.1572 

0.2084 

0.2400 

0.2549 

0.2553 

0,2621 

0.2990 

0.4407 

l.OOOO 

0.40 

-0. 

0.1440 

0.2773 

0.3935 

0-4878 

0-5561 

0.5925 

0,6071 

0.6085 

0-6493 

l.OOOO 

0.50 

-0. 

0. 1429 

0-2865 

0.4244 

0.5487 

0.6518 

0.7234 

0.7651 

0.7735 

0.7806 

l.OOOO 

0.60 

-0. 

0. 1221 

0.2555 

0.3940 

0.5292 

0-6519 

0.7485 

0,8144 

0,8383 

0.8386 

1,0000 

l.io 

-0. 

0.0416 

0.1128 

0.2117 

0.3352 

0.4773 

0.6240 

0.7628 

0.8845 

0-9663 

1.0000 

2.00 

-0. 

0.0421 

0.0879 

0.1573 

0.2566 

0.3796 

0.5235 

0.6790 

0.8249 

0.9421 

1.0000 

2.40 

-0. 

0.0472 

0-0920 

0.1554 

0-2481 

0.3654 

0.5040 

0-6568 

0,8065 

0.9283 

l.OOOO 

2.80 

-* o . 

0.0521 

0.0956 

0.1547 

0.2417 

0.3532 

0.4900 

0-6449 

0.7950 

0,9187 

l.OOOO 

3.20 

-0. 

0.0557 

0. 0994 

0.156 1 

0.2395 

0.3475 

0.4809 

0.6327 

0.7824 

0,9097 

1.0000 

4.00 

-0. 

0.0653 

0. 1135 

0.1660 

0,2418 

0.3426 

0.4696 

0-6174 

0.7646 

0,8911 

l.OOOO 

5.00 

-0. 

0.0936 

0. 1478 

0- 1968 

0.2660 

0.3578 

0.4756 

0.6149 

0.7533 

0,8741 

l.OOOO 

6.70 

-0. 

0. 1416 

0-2007 

0.243 1 

0,3003 

0.3787 

0.4780 

0.5952 

0.7136 

0.8237 

l.OOOO 

8.00 

~0. 

0. 1638 

0-2261 

0.2656 

0,3164 

0.3823 

0-4651 

0-5631 

0.6623 

0-7603 

l.OOOO 

9.00 

-0. 

0.1596 

0.2273 

0.2677 

0.3142 

0.3738 

0.4471 

0.5320 

0-6249 

0.7343 

1,0000 

10.30 

-0. 

0.1523 

0.2155 

0-2531 

0.2977 

0.3533 

0.4163 

0.4836 

0-5597 

0-6852 

1,0000 

12.50 

-0. 

0-1541 

0.2178 

0.2570 

0.2958 

0.3397 

0.3892 

0.4454 

0-5186 

0.6433 

1,0000 


♦01* EXIT IN ARGON 


TABLE XI. - INPUT TO PROGRAM ARGINV 


l 


cos 0 


Q -( cos "^ t , E)dt 
(t{E) 


SDATA 



- 1.0 

- 0.8 

- 0.6 

“ 0.4 

- 0.2 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

0.01 

0.0 

.0816 

.1661 

.2536 

.3446 

.4392 

.5378 

.6417 

.7517 

.8693 

1.0 

0.05 

0.0 

.0520 

.1098 

.1745 

.2471 

.3290 

. 4219 . 

.5289 

.6538 

.8033 

1.0 

0.10 

0.0 

.0235 

.0535 

.0918 

.1406 

.2029 

.2824 

.3862 

.5223 

.7081 

1.0 

0.20 

0.0 

.0041 

.0052 

.0053 

.0079 

.0192 

.0466 

.1094 

.2289 

.4601 

1.0 

0.30 

0.0 

.0870 

.1572 

.2084 

.2400 

.2 549 

.2553 

.2621 

.2990 

.4407 

1.0 

0.40 

0.0 

.1440 

.2773 

.3935 

.4878 

.5561 

.5925 

.6071 

.6085 

.6493 

1.0 

0.50 

c.o 

.1429 

.2865 

.4244 

.5487 

.6518 

.7234 

.7651 

.7735 

.7806 

1.0 

0.60 

0.0 

.1221 

.2555 

.3940 

.5292 

.6519 

.7485 

.8144 

.8383 

.8386 

1.0 

1.10 

0.0 

.0416 

.1128 

.2117 

.3352 

.4773 

.6240 

.7628 

.8845 

.9663 

1.0 

2.00 

0.0 

.0421 

.0879 

.1573 

.2566 

.3796 

.5235 

.6790 

.8249 

.9421 

1.0 

2.40 

0.0 

.0472 

.0920 

.1554 

.2481 

.3654 

.5040 

.6568 

.8065 

.9283 

1.0 

2.80 

0.0 

.0521 

.0956 

.1547 

.2417 

.3532 

.4900 

.6449 

.7950 

.9107 

1.0 

3.20 

0.0 

.0557 

.0994 

.1561 

.2395 

.3475 

.4809 

.6327 

.7824 

.9097 

1.0 

4.00 

o.c 

.0653 

.1135 

.1660 

.2418 

.3426 

.4696 

.6174 

.7646 

.8911 

1.0 

5.00 

c.o 

.0936 

.1478 

.1968 

.2660 

.3578 

.4756 

.6149 

.7533 

.8741 

1.0 

6.70 

0.0 

.1416 

.2007 

.2431 

.3003 

.3787 

.4780 

.5952 

.7136 

.8237 

1.0 

8.00 

0.0 

.1638 

.2261 

.2656 

.3164 

.3823 

.4651 

.5631 

• 6623 

.7603 

1.0 

9.00 

0.0 

.1596 

.2273 

.2677 

.3142 

.3738 

.4471 

.5320 

.6249 

.7343 

1.0 

10.3 

0.0 

• 1523 

.2155 

.2531 

.2977 

.3533 

.4163 

.4836 

.5597 

• 6852 

1.0 

12.5 

0.0 

.1541 

.2178 

.2570 

.2958 

.3397 

.3892 

.4454 

• 5186 

.6433 

1.0 


SIGMA-S 

3.713 
1.791 
0.901 
0.303 
0.197 
0.226 
0.276 
0.308 
0.270 
0.538 
0.656 
0.762 
0.8 8A 
1.116 
1.447 
2.053 
2.399 
7.849 
3.330 
3.531 
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TABLE Xn. - FIRST THREE PAGES OF PRINTED OUTPUT FROM ARCHNV (R,EV,COS 0) 


I R 


1 

0.0078125 

0.0234375 

0.0390625 

0.0546875 

0.0703125 

0.0859375 

0. 1015625 

0.1171875 

9 

0.1328125 

0.1484375 

0.1640625 

0.1796875 

0.1953125 

0.2109375 

0.2265625 

0.2421875 

17 

0.2578125 

0.27343 75 

0.2890625 

0.3046875 

0.3203125 

0.3359375 

0.3515625 

0.3671875 

25 

0.3828125 

0. 3984375 

0.4140625 

0.42968 75 

0.4453125 

0.4609375 

0.4765625 

0.4921875 

33 

0.5078125 

0.5234375 

0.5390625 

0.5546875 

0.5703125 

0.5859375 

0.6015625 

0.6171875 

41 

0.6328125 

0.6484375 

0.6640625 

0.6 796 875 

0.6953125 

0.7109375 

0.7265625 

0.7421875 

49 

0. 7578125 

0.7734375 

0.7890625 

0.8046875 

0.8203125 

0.8359375 

0.8515625 

0.8671875 

57 

0.8828125 

0.8984375 

0.9140625 

0.9296875 

0.9453125 

0.9609375 

0.9765625 

0.9921875 


J 

1 

0.0312500 

0.0937500 

0. 1562500 

0.2187500 

EV 

0.2812500 

0.3437500 

0.4062500 

0.4687500 

9 

0. 5312500 

0.5937500 

0.6562 5 00 

0.7187500 

0.7812500 

0.8437500 

0.9062500 

0.9687500 

17 

1.0312500 

1.0937500 

1.1562500 

1.2187500 

1.3750000 

1.6250000 

1.8750000 

2.1250000 

25 

2.3750000 

2.6250000 

2. 8750000 

3. 1250000 

3.3750000 

3.6250000 

3.8750000 

4.1250000 

33 

4.3750000 

4.6250000 

4. 8750000 

5.1250000 

5.3750000 

5.6250000 

5.8750000 

6.1250000 

41 

6.3750000 

6.6250000 

6.8750000 

7.1250000 

7.3750000 

7.6250000 

7.8750000 

8.1250000 

49 

8. 3750000 

8.6250000 

8. 8750000 

9.1250000 

9.3750000 

9.6250000 

9.8750000 

10.1250000 

57 

10.3750000 

10.6250000 

10.8750000 

11.1250000 

11.3750000 

11.6250000 

11.8750000 

12.1250000 


I J 


1 

1 

-0.9757 

-0.9387 

-0.5294 

-0.8949 

-0.9785 

-0.9873 

-0.9894 

-0.9894 

1 

9 

-0.9885 

-0.9871 

-0. 9853 

-0.9331 

-0.9805 

-0.9772 

-0.9733 

-0.9685 

1 

17 

-0.9630 

-0.9567 

-0.9503 

-0.9443 

-0.9358 

-0.9450 

-0.9589 

-0.9654 

1 

25 

-0.96 77 

-0.9701 

-0. 9721 

-0.9734 

-0.9744 

-0.9754 

-0.9768 

-0.9786 

1 

33 

-0.9806 

-0.9825 

-0.9842 

-0.9856 

-0.9867 

-0.9877 

-0.9885 

-0.9891 

1 

41 

-0.9897 

-0.9902 

-0.9907 

-0.9910 

-0.9913 

-0.9916 

-0.9917 

-0.9917 

1 

49 

-0.9917 

-0.9917 

-0.9915 

-0.9914 

-0.9913 

-0.9913 

-0,9911 

-0.9911 

1 

57 

-0.9911 

-0.9910 

-0.9910 

-0.9910 

-0.9910 

-0.9910 

-0.9910 

-0.9910 


2 

1 

-0.9273 

-0.8225 

-0.2829 

-0.5185 

-0.935 1 

-0.9620 

-0.9681 

-0.9681 

2 

9 

-0.9655 

-0.9613 

-0.9559 

-0.9495 

-0.9417 

-0.9323 

-0.9210 

-0.9079 

2 

17 

-0. 8930 

-0.8773 

-0.8620 

-0,8486 

-0,8291 

-0.8432 

-0.8774 

-0.8958 

2 

25 

-0.9025 

-0.9095 • 

-0.9156 

-0.9194 

-0.9224 

-0.9256 

-0.9298 

-0.9353 

2 

33 

-0.9413 

-0,9471 

-0.9523 

-0.95 65 

-0.9600 

-0.9628 

-0,9653 

-0.9673 

2 

41 

-0.9691 

-0.9706 

-0.9720 

-0.9731 

-0.9740 

-0.9747 

-0.9751 

-0.9752 

2 

49 

-0.9751 

-0.9748 

-0.9744 

-0.9740 

-0,9737 

-0.9735 

-0.9733 

-0.9731 

2 

57 

-0.9730 

-0.9730 

-0.9730 

-0.9730 

-0.9730 

-0.9730 

-0.9731 

-0.9731 


3 

1 

-0.8793 

-0.7200 

-0.1459 

0.0906 

-0.8909 

-0.9364 

-0.9469 

-0.9470 

3 

9 

-0.9426 

-0.9356 

-0.9268 

-0.9163 

-0.9037 

-0.8890 

-0.8718 

-0.8526 

3 

17 

-0.8320 

-0.8115 

-0. 7926 

-0.7764 

-0. 7509 

-0.7574 

“0.7976 

-0.8247 

3 

25 

-0.8353 

-0.8465 

-0. 8566 

-0.8631 

-0.8683 

-0.8740 

-0.8813 

-0.8907 

3 

33 

-0.9011 

-0.9110 

-0.9197 

-0.9269 

-0.9328 

-0.9376 

-0.9417 

-0.9452 

3 

41 

-0.9482 

-0.9508 

-0.9530 

-0.9550 

-0.9565 

-0,9576 

-0.9583 

-0.9585 

3 

49 

-0.9583 

-0.9578 

-0.9572 

-0.9566 

-0.9560 

-0.9556 

-0.9553 

-0.9550 

3 

57 

-0.9548 

-0.9547 

-0. 9547 

-0.9547 

-0.9547 

-0.9549 

-0.9550 

-0.9552 


4 

1 

-0.8319 

-0.6290 

-0.0470 

0.2212 

-0,8452 

-0.9108 

-0.9256 

-0.9258 

4 

9 

-0.9197 

-0.9100 

-0.8979 

-0.8837 

-0.8671 

-0.8480 

-0. 8266 

-0.8036 

4 

17 

-0.7801 

-0.7574 

-0. 7366 

-0.7137 

-0.6886 

-0.6861 

-0.7216 

-0.7517 

4 

25 

-0. 7644 

-0. 7786 

-0.7923 

-0.8018 

-0.8099 

-0.8187 

-0.8297 

-0.8435 

4 

33 

-0. 8586 

-0.8731 

-0. 8857 

-0.8963 

-0.9047 

-0.9117 

-0.9176 

-0.9226 

4 

41 

-0.9269 

-0.9306 

-0.9338 

-0.9365 

-0.9387 

-0.9403 

-0.9413 

-0.9416 

4 

49 

-0.9413 

-0.9406 

-0. 9397 

-0.9388 

-0.9381 

-0.9375 

-0.9370 

-0.9366 

4 

57 

-0.9364 

-0.9362 

-0.9362 

-0,93 62 

-0.9362 

-0.9365 

-0.9367 

-0.9369 
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TABLE Xn. - Concluded. FIRST THREE PAGES OF PRINTED OUTPUT FROM ARCHNV (R, EV, COS B) 

I J 


5 

1 

- 0.7352 

- 0.5464 

0.0324 

0.2985 

- 0.7969 

- 0.8849 

- 0.9042 

- 0.9046 

5 

9 

- 0.8968 

- 0.8845 

- 0. 8694 

- 0.8520 

- 0.8320 

- 0.8096 

- 0.7854 

- C .7602 

5 

17 

- 0.7349 

- 0.7105 

- 0.6883 

- 0.6690 

- 0.6353 

- 0.6257 

- 0.6534 

- 0.6804 

5 

25 

- 0.6919 

- 0.7054 

- 0. 7204 

- 0.7322 

- 0.7436 

- 0.7566 

- 0.7725 

- 0.7918 

5 

33 

- 0.8127 

- 0.8325 

- 0.8497 

- 0.8638 

- 0.8753 

- 0.8847 

- 0.8926 

- C .8992 

5 

41 

- 0.9049 

- 0.9098 

- 0.9140 

- 0.9176 

- 0.9205 

- 0.9226 

- 0.9239 

- 0.9243 

5 

49 

- 0.9239 

- 0.9231 

- 0.9219 

- 0.9208 

- 0.9198 

- 0.9189 

- 0.9183 

- 0.9178 

5 

57 

- 0.9174 

- 0.9172 

- 0.9172 

- 0.91 72 

- 0.9174 

- 0.9176 

- 0.9178 

- 0.9181 


6 

1 

- 0. 7394 

- 0.4705 

0. 0990 

0.3542 

- 0.7451 

- 0.3587 

- 0.8828 

- 0.8834 

6 

9 

- 0.8740 

- 0.8593 

- 0.8416 

- 0.8213 

- 0.7986 

- 0.7739 

- 0.7477 

- 0.7208 

6 

17 

- 0.6940 

- 0.6683 

- 0.6448 

- 0.6243 

- 0.5879 

- 0,5732 

- 0.5942 

- 0.6158 

6 

25 

- 0.6238 

- 0.6337 

- 0.6462 

- 0,6578 

- 0.6710 

- 0.6876 

- 0.7082 

- 0.7332 

6 

33 

- 0.7607 

- 0.7872 

- 0.8101 

- 0.8288 

- 0.8438 

- 0.8560 

- 0.8661 

- 0.8747 

6 

41 

- 0. 8820 

- 0.8882 

- 0.8936 

- 0.8981 

- 0.9017 

- 0.9044 

- 0.9060 

- 0.9065 

6 

49 

- 0.9061 

- 0.9050 

- 0.9036 

- 0.9021 

- 0.9009 

- 0.8998 

- 0.8990 

- 0.8984 

6 

57 

- 0.8979 

- 0.8976 

- 0.8976 

- 0.8976 

- 0.8978 

- 0.8980 

- 0.8984 

- 0.8988 


7 

1 

- 0,6945 

- 0.4005 

0.1560 

0. 3984 

- 0.6886 

- 0.8321 

- 0 . R 611 

- 0.8622 

7 

9 

- 0.8513 

- 0.8344 

- 0. 8143 

- 0.79 17 

- 0.7669 

- 0.7404 

- 0. 7126 

- 0.6843 

7 

17 

- 0.6561 

- 0.6292 

- 0.6047 

- 0.5832 

- 0.5448 

- 0.5266 

- 0.5426 

- 0.5592 

7 

25 

- 0.5633 

- 0.5692 

- 0.5780 

- 0.5873 

- 0.5996 

- 0.6165 

- 0.6391 

- 0.6678 

7 

33 

- 0, 7008 

- n . 7342 

- 0. 7643 

- 0.7891 

- 0.8088 

- 0.8247 

- 0.8376 

- 0.8484 

7 

41 

- 0.8575 

- 0.8654 

- 0. 8721 

- 0.8777 

- 0.8821 

- 0.8854 

- 0.8874 

- 0.8881 

7 

49 

- 0.8876 

- 0.8863 

- 0.8845 

- 0.8828 

- 0.8813 

- 0.8800 

- 0.8789 

- 0.8781 

7 

57 

- 0.8775 

- 0.8772 

- 0. 8770 

- 0.6770 

- 0.8773 

- 0.8776 

- 0.8781 

- 0.8786 


8 

1 

- 0.6502 

- 0.3356 

0.2059 

0.4356 

- 0.6264 

- 0.8050 

- 0.8395 

- 0.0409 

8 

9 

- 0. 8286 

- 0.8098 

- 0. 7876 

- 0.7631 

- 0.7366 

- 0.7086 

- 0.6794 

- 0.6498 

8 

17 

- 0.6205 

- 0.5926 

- 0.5671 

- 0.5449 

- 0.5050 

- 0.4843 

- 0.4965 

- 0.5091 

8 

25 

- 0. 5099 

- 0.5125 

- 0. 5183 

- 0.5252 

- 0.5350 

- 0.5498 

- 0.5707 

- 0.5992 

8 

33 

- 0.6343 

- 0.6728 

- 0. 7098 

- 0.7419 

- 0.7679 

- 0.7889 

- 0.8057 

- 0.8195 

8 

41 

- 0.8310 

- 0.8407 

- 0. 8490 

- 0.8560 

- 0.8614 

- 0.8654 

- 0.8679 

- 0.8687 

8 

49 

- 0.8681 

- 0. 8666 

- 0. 8645 

- 0.8625 

- 0.8606 

- 0.8590 

- 0.8577 

- 0.8567 

8 

57 

- 0.8559 

- 0.8555 

- 0.8553 

- 0.8553 

- 0.8556 

- 0.3560 

- 0.8566 

- 0.8572 
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TABLE Xm. - PRINTED OUTPUT FROM ENEC 


$PMl 


NO 

500, 

Nl 

= 

17, Kl = 

10, ALPHA = 

2.0000000E 00, CONST = 

5-OOOOOOOF 01 

NFLAG = 

2. 

BC 


-1.3000000E 01, KODE = 

?, MODE = 

1, ERROR = 

7. 9999999E-01 


NS = 

R, KFLAG = 

2, TEMPK = 

l.OOOOOOOE 03, LFLAG = 2, 


S ENC 

THFRMIONIC EMISSION 

152 




AI = 

-0.2I969125E 01 
O.6RO98A35E-02 
0. 13841715E-04 

0. 

-0.23765674F 00 
-0.208008AIE-02 
0.95622907F-04 

0. 

0.62491 354E 00 
0. 1 31143736-02 
-0.62333257F-04 

-0.15334123E 00 0. 5091 029 IF-O 1 

-0.829917236-03 0 . 57 1915 79F - 03 

0. 194052106-03 -0. 1 26632466- 04 

-0.199465246-01 

-0.43561954E-03 

0. 

XD 

MEAN 

N 

STD. N 



1 

-0. 

96 7,90 0 

4,537 




2 

0. 00960736 

856.800 

5.957 




3 

0.03806023 

613.200 

10.594 




4 

0.08426519 

390.200 

8.821 




5 

0.14644661 

253.600 

8. 384 




6 

0.22221488 

163.000 

7.467 




7 

0.308o5828 

1 1 l.OOO 

4.869 




B 

0.40245484 

76.600 

4.387 




9 

0-49999999 

61. BOO 

4.320 




10 

0. 59754515 

51.400 

4.248 




11 

0. 69134171 

45,000 

4.768 




12 

0. 777785U 

41.400 

5.012 




13 

0.85355338 

38.400 

5.431 




14 

0.91573481 

35.200 

4.791 




15 

0.96193976 

35.000 

4.879 




16 

C. 99039263 

33.000 

5.013 




17 

1. 00000000 

32.000 

4.556 





MEAN AI 

STD. AI 

MEAN DA 

STD, DA 

MEAN B 

STD. a 


1 

-3.66167679 

0.17469033 

-6.51143974 

0.39980942 

41.15168715 

0.60209867 

2 

- 1.04503354 

0.11007060 

5.43993092 

0, 16972015 

-30.33697963 

0.23735175 

3 

0.54961915 

0,02118890 

-2.33130553 

0.06242871 

19.39196372 

0.27361134 

4 

-0.15537757 

0.00624523 

1 .04297768 

0.04069808 

-11.6865351? 

0.39588887 

5 

0.05116830 

0.00363047 

-0.46677465 

0.02416211 

6.87623149 

0.48478547 

6 

-0.01727193 

0.0019508? 

0,22428477 

0.01B97346 

-4.21814072 

C. 36597438 

7 

0.00658920 

0,00065380 

-0.17133606 

0.018791C6 

2.39053601 

0. 19135751 

8 

-0.00266536 

0. 00080612 

0.06614400 

0.00935168 

-1.30607525 

0.24036396 

9 

0.00110793 

0.0005041 1 

-0,04670603 

0.005783C2 

0.53050392 

0-21 118020 

10 

-0.00122989 

0.00023203 

0.03069014 

0.00844622 

0. 1885178B 

0.25728086 

11 

0.00094552 

0.00036242 

-0.00242998 

0.00321652 

-0.56634136 

0. 32823963 

12 

-0,00014676 

0.00017450 

-0.00713085 

0.00795397 

0.28571755 

0.1933983? 

13 

0.00003092 

0. 000166 70 

0.00402765 

0.00709825 

-0.25258386 

0,20521291 

14 

-0.00016668 

0.000 16971 

-0.00861512 

0.00366031 

0.09239006 

0.24925443 

15 

-0.00005572 

0.00012669 

0,01269518 

0.00407287 

0,1954028? 

0.29722011 

16 

0.00037266 

0.00013171 

-0,00549471 

0.00404726 

-0.61854020 

0.30506335 

17 

-0.00008585 

0. 00006323 

-0.00966468 

0.00476661 

0,52508580 

0.13166420 
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TABLE xm. - Concluded. PRINTED OUTPUT FROM ENEC 



MEAN Y STO. Y 

MEAN DY 

STD. OY 

MEAN DDY 

STD. OOY 

1 

O.UOOOOOlO 0.00000001 

-13.00000048 

0.00000004 

97.27428341 

1-95648837 

2 

-0.12043346 0-00011S87 

-12.14068496 

0.01643634 

85,69345951 

1.59331208 

3 

-0.43435966 0-00082861 

-10,01470077 

0.03406554 

63,62157488 

1.25241616 

A 

-0.84042504 0,00225005 

-7.77618700 

0.05459954 

36.70868532 

0- 94069674 

5 

-1.26187769 0.00642534 

-5.87558651 

0-09958778 

26.71593428 

1.09402619 

6 

-1.63894701 0.01630090 

-4. 19454759 

0.18267000 

17-47496819 

1.12554373 

7 

-1.94463196 0.03347100 

-2.95429876 

0.21243631 

12.34195209 

0.76109815 

S 

-2.17403677 0,05201884 

-2.00364363 

0.20529954 

7-75743997 

0.51680201 

9 

-2.33571520 0-07211596 

-1. 33246121 

0-24044415 

6.85151982 

0-65472865 

10 

-2.43513528 0.09561784 

-0.73212643 

0. 27012233 

4.69202816 

0.63100371 

11 

-2.48597673 0.12031564 

-0.37907971 

0.28503805 

3.66198051 

0.34921714 

12 

-2.50515202 0.14432648 

-0.05743694 

0,29545915 

3.31606522 

0-31601627 

13 

-2.50015339 0.16649128 

0. 17689936 

0.31401107 

3-31563717 

0.67820865 

14 

-2.48312446 0.18585514 

0.37558117 

0.32977037 

2-45997602 

0.31569276 

15 

-2.46316043 0.20063088 

0.47336564 

0.33487796 

2.62418780 

0.39094478 

16 

-2.44876626 0.20988189 

0.56060471 

0-34028953 

2.22343302 

0.26049839 

17 

-2.44303799 0.21302953 

0.56557281 

0.34145060 

2.07499668 

0.22437462 


ME AN 

STD. 





CURRENT 0.06400000 

0.00911165 





VOLT -2.44303799 

0. 21302953 





OPHIO -13.00000048 

0. 00000004 





KNTR 702.100 

15.755 





PHIMIN= -2.50568336 






XMIN= 0. 794921 88 





EQUALLY 

SPACED DATA 






X 

PHM X» 

DPHItX) 

NtX ) 

0- 

0 . nononoi 3 

-13.00000036 

1.94548561 

0.05000000 

-0, 54952634 

-9. 31162477 

1.10364881 

0. 099 >99J9 

-0.95850401 

-7.22976363 

0.64567125 

0. 15000000 

-1.2B257124 

-5.78223091 

0.52781617 

'J. 2000UOOO 

-1.54111073 

-4,61686987 

0.40628800 

0. 25000000 

-1.74916059 

-3.74006915 

0.29716718 

0. 30000000 

-1.91863452 

-3.06029117 

0.25245698 

0. 34999999 

-2.05670783 

-2.48759082 

0.21 1 38898 

0.40000000 

-2. 16908213 

-2.02345943 

0 ,15734602 

0.45000030 

-2.26104784 

-1.66022669 

0.13258342 

0. 50000000 

-2. 33571520 

-1.33246116 

0.13703039 

0. 55000000 

-2. 39395094 

-1.00442715 

0.12625750 

0. 59999999 

-2.43692946 

-0. 72007295 

0.092 16407 

0.65000000 

-2. 46764050 

-0.5 1697391 

0.07136576 

0. 70000000 

-2. 489091 01 

-0. 34924910 

0.07421897 

0. 75000000 

-2.50197527 

-0.16180167 

0,072 1677? 

0. BOOUOOOO 

-2.50563982 

0.01771381 

0,06333854 

0. 84999999 

-2.50078127 

0, 16612880 

0.06619857 

0. 90000000 

-2.48856401 

0.32855473 

0.05539818 

0.95000000 

-2. 468764 04 

0.44869374 

0.04913158 

1,00000000 

-2.44303799 

0.56557281 

0.04149993 
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I 


X*F-0 Y*e 1 



0. O.lOO 0.200 0.300 

l.OOl 

0.400 0.500 0,600 

0.700 

0. ROO 

0.900 

1.000 


1 

1 

1 

1 







1 

1 






0.000 0. 

1 

0. ♦ 
1 






-0.550 0.500 

1 

1 

1 ♦ 

1 

1 






-0.959 1.000 

1 

1 

-l.OOl ♦ 

1 






-1.7R3 1.500 

1 

1 ♦ 






-1.541 7.000 

1 ♦ 






-1.749 2.500 

L * 






-1.919 3.000 

-2.057 3.500 
1A9 4.000 
-2.336 5.000 
-2.44 310.000 
-7.506 8.000 

1 ♦ 

-2.001 

1 * 
I 
1 
1 
1 
1 
1 

* 

* * 

* * 

* 

♦ ♦ 

* * 

* * 

• 


1 

1 







9. O.lOO 0.200 0.300 

0.400 0.500 0.600 

0. 700 

0. 800 

0.900 

l.OOO 


SCALE FACTOR FOR OROINATES IS 10**( 0) 

(a) Printer plot of ^(X). 





X*F-0 Y*F 1 

0, O.lOO 0.200 0,300 

1 00 1 - — - 

0.400 0,500 0.600 

C. 700 

0.800 

0.900 

l.OOO 


1 

1 

I 

1 

1 

1 

1 

1 

I 

0. 381 
1 
1 






0.033 9.000 
-0.016 7.500 
-0.07? 6.000 
-0. 133 5.000 
-0.20? 4.000 
-0.749 3.500 
-0.306 3.000 
-0.374 2.500 
-0.462 2.000 

1 

1 

1 

1 

1 

1 

1 

-0.251 ♦ 

1 ♦ 

1 ♦ 

1 + 

♦ ♦ 

* * 

* * 

* 

* 

* * 

* * 

♦ 

-0,578 1.500 

1 ♦ 
1 






-0.723 l.OOO 
-0.931 0.500 

1 

1 * 
1 

-0.871 

1 ♦ 

1 

1 







1 

1 






-1.300 0. 

♦ 







1 

-1 501 - - - - 







0. O.lOO 0.200 0.300 

0.400 0.500 0.600 

0,700 

0.800 

0.900 

1.000 


SCALE FACrriR FOR ORDINATES IS 10**( U 

(b) Printer plot of <o'(X). 
Figure 39. - Output from ENEC. 
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X*E~1 Y*E 1 


0.973 0. 


0.552 0.500 


0 . 0.100 

10.01 

* 

1 

1 

1 

1 

1 

I 

1 

I 

8.01 

1 

1 

1 

1 

1 

1 

1 

1 

I 

6.01 

I 

1 ♦ 


0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900 


I 

I 

I 

1 

4.01 

I 


0.323 1.000 


0.264 1.500 


0.203 2.000 


0.149 2.500 
0.126 3.000 
0.106 3.500 
0.079 4.000 
0.063 5.500 
0.032 8.000 
0.02110.000 


1 

1 

1 

1 

1 

1 

1 

2,01 

1 

1 

1 

1 

1 

1 

1 

1 

1 

0 . 1 

0 . 




♦ ♦ ♦ 


♦ ♦ ♦ ♦ 


0.100 0.200 0.300 


0.400 0.500 


0.600 


0,700 


0.800 0.900 


SCALE FACTOR FOR ORDINATES IS 104*1 1! 


(c) Printer plot of <p”(X). 


Figure 39. ~ Concluded. 
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*'The aeronautical and space activities of the United States shall be 
conducted so as to contribute ,,, to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof/^ 

— ^National Aeronautics and Space Act of 1958 
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